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EXECUTIVE  SUMMARY 


This  report  contains  one  interim  report,  together  with  the  most  recent  results.  Part  I  describes  the 
work  done  during  a  visit  by  Mr  Bosch  of  the  Topographic  Engineering  Center,  Alexandria.  Part  II  describes 
the  analysis  of  selected  pixels  from  part  of  a  SPOT  image  of  Fort  A.  P.  Hill.  They  were  chosen  to  represent 
several  different  types  of  ground  cover.  The  NIR  part  of  the  spectrum  was  analysed  for  these  pixels.  There 
were  seven  classes  of  pixels  and  the  values  in  each  were  subject  to  an  exploratory  data  analysis.  The  results 
show  that  several  classes  have  distinct  characteristics  that  were  distinguished  by  the  variograms  of  the  class 
data  and  their  model  parameters. 

Part  m  describes  the  analysis  of  hyperspectral  data:  the  hydice  imagery  for  Fort  Hood  and  hymap  imagery 
for  Fort  A.  P.  Hill.  Pixels  were  selected  for  six  different  ground  cover  classes  at  each  site.  The  spectrum  for 
each  pixel  was  plotted  and  examined  visually.  All  of  the  experimental  variograms  computed  for  each 
spectrum  showed  evidence  of  trend  and  this  was  usually  removed  satisfactorily  by  linear  and  quadratic 
functions.  The  variograms  were  recomputed  on  the  residuals  and  modelled  with  several  periodic  functions 
for  which  a  program  was  written  in  GenStat.  For  Fort  Hood  the  spectra  vary  considerable  for  the  different 
ground  cover  classes  and  the  variograms  have  certain  distinctive  features.  The  spectra  for  Fort  A.  P.  Hill  are 
less  distinctive  apart  from  one  class,  but  the  amplitudes  of  the  periodic  functions  fitted  to  the  variograms 
provide  one  means  of  discrimination.  To  explore  these  data  further  principal  component  analyses  were 
done.  For  both  sites  three  components  accounted  for  most  of  the  variation.  Certain  sections  of  the  spectrum 
accounted  for  the  variation  on  these  axes.  When  the  pixels  were  plotted  in  a  projection  of  the  first  two 
components  four  clear  groups  were  evident  for  Fort  Hood  and  three  for  A.  P.  Hill.  A  non-hierarchical 
multivariate  method  of  classification  was  applied  to  50%  of  the  waveband  data  and  to  the  scores  of  the 
leading  principal  components.  The  results  using  the  latter  provided  the  best  discrimination  and  the  best 
relation  with  the  ground  cover  classes.  There  are  features  in  the  spectra  and  the  variograms  that  enable 
discrimination  between  ground  cover  classes,  which  is  supported  by  the  multivariate  analyses. 

Appendix  I  is  a  GenStat  program  to  fit  periodic  functions  to  the  experimental  variogram.  In  addition 
Appendix  II  is  the  text  of  another  paper  resulting  from  work  with  Mr  Shine  and  Mr  Slocum.  This 
has  been  submitted  to  the  International  Journal  of  Remote  Sensing.  The  comparison  between  the 
usual  variogram  estimator  (Matheron’s  method  of  moments)  and  the  wavelet  based  variograms  is 
given  in  Appendix  DL  The  results  suggest  that  the  wavelet  variograms  can  remove  long-range  trend 
features  whether  they  arise  from  stepped  functions,  such  as  those  related  to  the  boundaries  between 
different  types  of  soil  or  rock,  or  from  smooth  variation  that  can  be  represented  by  a  trend  surface. 
This  would  avoid  the  need  to  examine  the  data  for  trend  first.  The  effects  of  local  trend  or  drift, 
which  are  less  common  in  most  earth  science  data,  appear  to  be  amplified  by  the  wavelet-based 
variograms,  however.  An  oral  presentation  of  this  work  will  be  made  at  the  GeostatsUK  meeting  in 
Reading,  England,  in  March  2002.  Professor  Webster  has  devoted  most  of  his  time  on  this  contract 
to  preparing  the  theoretical  part  of  the  text  with  Mr  Bosch  for  publication  and  writing  Fortran  code 
for  the  simulations. 


PART  I  REPORT  ON  VISIT  BY  E.  BOSCH 


The  purpose  of  the  visit  by  Mr  Bosch  was  to  continue  with  our  work  to  publish  a 
paper  on  the  wavelet  variogram.  This  was  started  as  part  of  the  previous  project  and  is 
now  nearing  completion.  Professor  Webster  has  made  considerable  input  into  this  in 
making  the  mathematics  more  comprehensible  to  non-mathematicians.  The  paper 
shows  the  relation  between  the  conventional  variogram  of  geostatistics  and  a 
variogram  based  on  the  wavelet  transform.  The  discovery  of  the  relation  by  Mr  Bosch 
is  novel  as  far  as  we  know. 

Professor  Webster  spent  three  days  working  with  us  on  this  and  during  this  time  many 
analyses  were  done  to  try  to  gain  insight  into  what  the  two  approaches  do.  The 
analyses  were  done  on  two  sets  of  real  data  and  one  set  of  simulated  data.  It  seems 
that  if  there  is  long-range  trend  present  in  the  data  the  wavelet  variogram  can  deal 
with  it  whereas  it  must  be  removed  before  the  analysis  of  the  conventional  variogram. 
Also  where  there  are  distinct  boundaries  (different  classes  of  soil  along  a  transect)  it 
seems  that  the  Daubechies  wavelet  variogram  disregards  the  effect  of  the  boundaries 
and  the  variogram  relates  to  the  variation  between  the  class  boundaries. 

The  other  two  days  that  Mr  Bosch  was  in  Reading  were  spent  in  preparing  data  for 
analysis  during  the  second  part  of  this  project.  In  addition,  we  examined  different 
levels  of  thresholding  in  a  wavelet  analysis  of  the  Fort  A.  P.  Hill  SPOT  data  The  aim 
of  this  is  to  compare  wavelet  multiresolution  analysis  and  factorial  kriging.  These 
results  will  be  added  to  the  paper  presented  at  the  International  Geostatistics  Congress 
in  2000  for  publication  in  the  International  Journal  of  Remote  Sensing. 


PART  II  ANALYSIS  OF  SELECTED  CLASSES  OF  GROUND  COVER  AT 

FORT  A.  P.  HILL 

Figure  la  shows  the  NIR  values  from  the  SPOT  image  of  Fort  A.  P.  Hill  as  a  pixel 
map.  Figure  lb,  derived  from  Matlab,  shows  the  location  of  the  seven  classes.  By 
comparing  the  two  maps  it  is  clear  that  the  classes  relate  to  areas  with  identifiable 
levels  of  reflectance  in  the  NIR  waveband.  Some  classes  are  similar,  but  they  are  in 
different  parts  of  the  scene.  Class  10  is  the  least  homogeneous,  which  can  be  seen  if 
Figure  la  is  examined  carefully.  It  incorporates  both  valley  and  spur  components. 
Class  20  appears  to  be  an  open  area  from  the  topographic  sheet  and  on  a  valley  side 
between  a  spur  or  upland  area  and  valley  bottom.  Classes  30  and  40  are  upland  or  spur 
areas  with  low  to  intermediate  reflectance  values.  Classes  50  and  60  are  valley  areas 
and  appear  to  be  more  homogeneous  than  many  of  the  other  classes.  Class  70  appears 
to  be  unwooded  in  part  and  on  a  spur. 

Table  1  gives  the  summary  statistics  of  the  seven  ground  cover  classes  for  near  infra 
red  (NIR)spectral  values.  Class  10  is  the  most  variable  class  (see  Figure  2a)  and  also 
in  terms  of  the  coefficient  of  variation.  Classes  20  and  40  are  the  least  variable,  i.e.  the 
most  homogeneous  (see  Figures  3a  and  5a).  The  remaining  four  classes  have  similar 
coefficients  of  variation.  Table  1.  The  less  variation  there  is  within  the  class  the 
greater  its  coherence.  Classes  20  and  30  have  similar  mean  values,  as  do  classes  40 
and  70,  and  50  and  60.  Class  10,  as  already  noted,  is  different  from  the  other  classes 
and  this  is  also  the  case  with  respect  to  the  mean. 
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The  histogram  of  class  10  suggests  that  the  statistical  distribution  is  bimodal.  Again 
this  shows  its  difference  from  the  other  classes  and  probably  reflects  the  spur  and 
valley  components  in  the  class.  As  a  consequence  it  might  have  been  beneficial  to  do 
the  variogram  analysis  of  the  two  strata  separately.  However,  there  were  insufficient 
data  for  this.  All  of  the  other  histograms  show  a  unimodal  distribution.  The  values  for 
all  classes  have  a  near-normal  distribution  as  shown  by  the  histograms  and  the 
skewness  values.  Therefore,  no  transformation  of  the  data  was  necessary. 


Table  1:  Summary  statistics  for  the  six  spectral  classes  for  Fort  A.  P.  Hill. 


Classes 


Parameter 

10 

20 

30 

40 

50 

60 

70 

Number  of 
observations 

183 

146 

192 

161 

167 

193 

113 

Mean 

126.0 

118.1 

114.0 

98.59 

141.1 

140.1 

102.2 

Median 

128.0 

119.0 

113.0 

99.0 

142.0 

141.0 

103.0 

Minimum 

107.0 

109.0 

95.0 

90.0 

119.0 

122.0 

90.0 

Maximum 

148.0 

125.0 

131.0 

110.0 

166.0 

169.0 

116.0 

Variance 

102.9 

13.62 

60.67 

15.33 

78.06 

68.03 

28.38 

Standard 

deviation 

10.14 

3.690 

7.789 

3.915 

8.835 

8.248 

5.327 

Coefficient 
of  variation 

8.048 

3.127 

6.832 

3.971 

6.266 

5.891 

5.223 

Skewness 

-0.160 

-0.443 

0.237 

0.115 

0.085 

0.359 

0.024 

Kurtosis 

-1.172 

-0.134 

-0.606 

-0.098 

0.053 

0.934 

-0.046 
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Spot  AP  Hill  Masked  Areas 


Figure  1.  a)  Part  of  SPOT  image  for  Fort  A.  P.  Hill,  b)  locations  of  six  of  the  seven 
classes  from  the  SPOT  image  of  A.  P.  Hill. 
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Class  1 0 

Figure  2a  shows  the  pixel  map  of  the  NIR  values,  which  cover  a  wide  range.  There  is 
an  upland  area  in  the  middle  between  two  lowland  or  valley  areas.  The  axes  are 
pixels,  each  of  which  is  20  m  x  20  m.  The  variogram  for  class  10  was  the  most 
difficult  to  model.  It  is  upwardly  concave  near  to  the  origin  suggesting  the  presence  of 
local  drift.  A  quadratic  function  fitted  on  the  coordinates  accounted  for  46%  of  the 
variation.  However,  when  the  variogram  was  computed  afresh  on  the  residuals  from 
the  quadratic  trend  its  shape  remained  the  same.  Therefore,  the  analysis  was  done  on 
the  raw  data.  The  best  fitting  models  were  the  circular  and  the  stable  exponential 
function.  The  latter  aims  to  avoid  the  disadvantages  of  the  Gaussian  model.  Visually 
the  range  of  spatial  dependence  appears  to  be  a  little  less  than  199m  as  indicated  by 
the  circular  model.  The  range  for  the  stable  exponential  model  is  not  sensible. 
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Figure  2.  Class  10:  a)  pixel  values  of  NIR,  b)  histogram,  c)  experimental  variogram 
(symbols)  circular  variogram  model,  (line),  d)  experimental  variogram  (symbols) 
stable  exponential  variogram  model. 
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Class  20 

Figure  3a  shows  the  NIR  values  plotted  as  a  pixel  map  for  class  20.  It  seems  that  this 
class  might  comprise  a  clearing  between  a  valley  and  spur.  The  NIR  values  are  in  the 
upper  part  of  the  range.  As  before  the  pixels  are  20  m  x  20  m  Figure  3c  shows  the 
experimental  variogram  and  Figure  3d  the  fitted  circular  model.  It.  has  a  range  of 
spatial  dependence  similar  to  that  for  class  10.  The  sill  variance,  however,  is  much 
smaller  showing  the  effects  of  the  smaller  maximum  reflectance  values  and  smaller 
degree  of  variation  present. 


Class  20 


Circular  model 


Lag  distance  /  m 


Figure  3.  Class  20:  a)  pixel  values  of  NIR,  b)  histogram,  c)  experimental  variogram 
(symbols),  d)  experimental  variogram  (symbols)  circular  variogram  model  (line). 
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Class  30 

Figure  4a  of  the  pixel  values  shows  that  this  class  has  a  markedly  patchy  variation.  It 
is  largely  an  upland  area,  but  possibly  with  small  patches  or  lower  lying  areas  within 
it.  The  NIR  values  are  intermediate  to  large.  The  best  fitting  model  was  a  nested 
spherical  one  with  a  short-range  structure  of  about  68  m  and  a  long-range  structure  of 
about  175  m.  The  latter  is  similar  to  the  ranges  for  classes  10  and  20.  This  class  is 
distinct  from  all  of  the  others  in  having  this  nested  variation  which  is  typical  of  the 
whole  image,  but  with  larger  ranges. 


Class  30 


NIR 


Double  spherical  model 


Figure  4.  Class  30:  a)  pixel  values  of  NIR,  b)  histogram,  c)  experimental  variogram 
(symbols,  d)  experimental  variogram  (symbols)  double  spherical  variogram  model 
(line). 
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Class  40 

Class  40  appears  to  be  the  most  uniform  (Figure  5a),  however,  its  coefficient  of 
variation  is  larger  than  for  class  20.  The  NIR  values  are  in  the  small  to  intermediate 
range.  It  is  an  upland  area  that  appears  to  have  been  cleared  of  woodland.  The 
variogram  has  the  smallest  sill  variance,  which  supports  the  visual  impression  gained 
from  this  class.  The  variogram  was  fitted  best  by  a  pentaspherical  model.  Although 
the  variation  appeared  to  be  the  least,  the  variogram  of  this  class  has  the  shortest  range 
of  spatial  dependence  of  all  variograms.  Table  2. 


Class  40 


NIR 


Penta  spherical  model 


Figure  5.  Class  40:  a)  pixel  values  of  NIR,  b)  histogram,  c)  experimental  variogram 
(symbols),  d)  experimental  variogram  (symbols)  pentaspherical  variogram  model 
(line). 
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Class  50 

This  class  also  appears  similar  visually  (Figure  6a)  to  class  60,  although  it  is  slightly 
more  variable.  They  are  both  in  valley  locations.  The  NIR  values  are  in  the  upper 
range.  The  sill  variance  of  the  variogram  is  relatively  large  and  the  range  of  spatial 
dependence  is  different  from  all  other  classes.  It  is  intermediate  at  1 1 7  m. 


Class  50 


Figure  6.  Class  50:  a)  pixel  values  of  NIR,  b)  histogram,  c)  experimental  variogram 
(symbols),  d)  experimental  variogram  (symbols)  circular  variogram  model  (line). 


Class  60 

Figure  7a  shows  the  pixel  values  for  class  60,  which  are  in  a  similar  range  to  those  for 
class  50.  It  is  a  valley  location.  The  variogram  has  an  intermediate  sill  variance,  which 
is  similar  to  that  for  class  50.  Its  range  of  spatial  dependence  of  103  m  is  the  second 
smallest.  Therefore,  although  the  pixel  values  are  similar  to  class  50  the  texture  in  the 
image  defined  by  the  range  is  different. 


Class  60 


Circular  model 


*  x 


X 


x  X 


X 


X 


50 

40 


X 


30 

20 

10  H 


”1 - 1 - 1 - 1 - f“ 

50  75  100  125  150 


“I - T“ 

175  200 


Figure  7.  Class  60:  a)  pixel  values  of  NIR,  b)  histogram,  c)  experimental  variogram 
(symbols),  d)  experimental  variogram  (symbols)  circular  variogram  model  (line). 
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Class  70 

Figure  8a  shows  the  pixel  map  of  NIR  values,  which  are  small  to  intermediate  and 
similar  to  those  of  class  40.  This  class  is  in  an  upland  area  and  there  might  be  some 
unwooded  area  that  give  rise  to  the  evident  patchiness  in  the  spatial  distribution  of 
NIR  values.  The  variogram  sill  is  modest  and  of  a  similar  order  of  magnitude  to  those 
for  classes  20  and  40,  which  are  in  similar  landscape  positions.  The  variogram  range 
of  139  m  is  intermediate  and  different  from  all  of  the  other  classes. 


Class  70 


Figure  8.  Class  70:  a)  pixel  values  of  NIR,  b)  histogram,  c)  experimental  variogram 
(symbols),  d)  experimental  variogram  (symbols)  pentaspherical  variogram  model 
(line). 
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Table  2:  Model  parameters  for  the  functions  fitted  to  the  experimental  variograms  of 

the  six  spectral  classes  for  Fort  A.  P.  Hill. 


Model  parameter 


Class 

Model  type 

Nugget 

variance 

Sill  (1) 
variance 

Sill  (2) 

variance 

Range  1 
(m) 

Range  2 
(m) 

Exponent 

Class  10 

Circular 

0 

163.3 

199.0 

Stable 

exponential 

0 

165.1 

309.9 

1.973 

Class  20 

Circular 

0 

20.13 

180.9 

Class  30 

Double 

spherical 

0.1479 

32.05 

33225 

67.92 

175.6 

Class  40 

Penta- 

spherical 

3.288 

12.36 

90.86 

Class  50 

Circular 

0 

95.31 

117.4 

Class  60 

Circular 

0 

74.46 

102.6 

Class  70 

Penta- 

spherical 

3.714 

22.03 

139.4 

Summary  of  analysis  of  classes 

The  variogram  ranges  identified  for  these  classes  relate  to  ones  already  described  for 
the  SPOT  image  and  the  1-m  CAMIS  imagery.  The  short-range  component  of  the 
variation  identified  in  the  SPOT  image  was  130  m,  which  relates  to  the  variation 
described  by  class  70.  The  long-range  component  of  420  m  was  not  present  in  any 
class  which  is  as  expected.  The  long-range  component  in  the  SPOT  image  represents 
the  main  physiographic  elements  and  associated  differences  in  ground  cover.  The 
short-range  component  identified  in  the  variogram  of  class  30  of  about  68  m  was 
identified  in  the  NIR  waveband  of  the  1  in  15  selection  from  the  1-m  imagery  of  75  m 
for  the  average  variogram  (see  Table  12  report  N68171-97-C-9029).  The  group  of 
ranges  of  about  1 00  m  and  1 90  m  were  also  detected  in  the  1  -m  imagery. 

The  results  suggest  that  parts  of  the  image  associated  with  particular  types  of  ground 
cover  have  a  texture  that  can  be  described  by  the  variogram  ranges  and  the  sill 
variances.  Classes  30,  50  and  70  are  different  from  each  other  and  the  other  classes  in 
terms  of  their  sill  variances  and  ranges  of  the  variogram  models.  Classes  10  and  20 
have  similar  ranges  but  very  different  variances;  classes  40  and  60  have  similar 
ranges,  but  again  very  different  sill  variances,  reflectance  values  and  they  are  in 
different  landscape  positions.  Classes  50  and  60  have  fairly  similar  variograms  and 
reflectance  values,  but  the  variogram  ranges  are  different.  Classes  40  and  70  have 
similar  nuggetrsill  ratios,  similar  reflectance  values  and  landscape  attributes. 
Nevertheless  it  seems  that  each  of  these  classes  can  be  distinguished  from  other 
classes  either  in  terms  of  the  variogram  range  or  the  sill  variance. 
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PART  III  ANALYSIS  OF  HYPERSPECTRAL  IMAGERY  FROM  FORT 
HOOD  AND  FORT  A.  P.  HILL 

The  aims  of  this  investigation  were  to  examine  the  spectra  of  individual  pixels  from 
different  types  of  ground  cover  to  assess  their  similarity  or  otherwise,  and  to 
determine  whether  their  variograms  were  distinctive.  In  addition  a  principal 
component  analysis  was  done  together  with  non-hierarchical  classifications  using  the 
spectral  values  and  the  principal  component  scores.  Spectral  information  from  Fort 
Hood  and  Fort  A.  P.  Hill  were  examined. 

Fort  Hood  spectra 

Pixels  from  a  Hydice  hyperspectral  image  of  Fort  Hood  (Texas)  near  to  Clean  were 
selected  according  to  different  types  of  ground  cover  comprising  different  types  of 
vegetation,  asphalt  and  soil.  The  image  is  shown  in  Figure  9a.  The  3  m  x  3  m  pixels 
in  each  of  the  six  groups  were  contiguous  and  relate  to  the  areas  marked  by  the 
arrows.  Figure  9b  shows  the  relative  sizes  of  the  classes  of  pixels  and  their  locations. 
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Figure  9.  a)  Colour  image  of  part  of  Fort  Hood  and  b)  the  locations  and  sizes  of  the 
classes  of  pixels. 


There  were  132  wavelength  values  (nm)  in  each  spectrum  and  a  total  of  329  pixels. 
Table  3  gives  the  number  of  sites  in  each  class  and  the  associated  ground  cover  class. 
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Table  3:  Classes  of  different  ground  cover  type  and  number  of  pixels  in  each  for  Fort 
Hood. 


Class 

Type  of  ground 
cover 

Number  of 
pixels 

Pixel  range 

1 

Grass 

70 

1-  70 

2 

Trees 

92 

71  -  162 

3 

Asphalt 

66 

163-228 

4 

Bare  soil 

21 

229  -  249 

5 

Mixed  herbaceous 
vegetation 

38 

250  -  287 

6 

Disturbed  soil 

42 

288-329 

The  spectrum  for  each  pixel  was  plotted  as  a  trace  and  examined  for  differences  and 
similarities.  Experimental  variograms  were  computed  from  the  spectra  for  each  pixel. 
The  full  set  of  spectra  and  variograms  can  be  provided,  but  to  simplify  interpretation 
we  focus  on  selected  ones  with  particular  properties.  This  selection  was  based  on  the 
form  of  the  spectrum  and  its  variogram.  This  part  of  the  investigation  was  done 
without  reference  to  the  classes  of  ground  cover  in  case  there  were  differences  in  the 
spectra  for  a  given  classes  which  might  be  overlooked  by  focusing  on  the  classes 
only. 


Spectra 

The  selected  spectra  for  the  six  ground  cover  classes  are  shown  in  Figures  12  to  84. 
Some  spectra  have  distinctly  different  forms,  while  others  vary  in  the  finer  detail  or  in 
their  maximum  reflectance  values. 

A  visual  inspection  of  the  spectra  suggested  the  following  similarities  and  differences 
for  Fort  Hood.  The  lists  marked  by  G  are  spectra  belonging  to  a  pure  and  distinctive 
group,  those  marked  by  -fc  are  similar  although  they  belong  to  different  ground  cover 
classes. 
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Pixels: 


G  1-70  (Class  1) 

71-82 

f  83  -  84,  86  -  89,  94,  96  -  98,  103,  105  -  107,  1 12  -  1 13,  1 15  -  1 19,  121  -  123, 


125 


-130, 

146-  152,  156,  158-  160,  162. 
250  -  264,  266  -  271,  273  -  283, 


285  -  287,  265,  272,  284.  (Class  5) 


85,  93,  102. 

90,  111. 

91  -92,  99-  101 
95 

104,  108  -  110,  114, 120,  124,  143  -  145,  153  -  155,  157,  161. 


>  (Class  2) 


131  -  133 
134-139 

140 

141  -  142 

G  163  -228  (Class  3) 

G  229  -  249  (Class  4) 

G  288  -  325,  328  -  329  (Class  6) 
326  -  327 


J 


Variography 

All  of  the  experimental  variograms  showed  evidence  of  trend  with  an  upwardly 
concave  slope  near  to  the  origin.  Linear  and  quadratic  functions  were  fitted  to  the 
wavelengths  (one  dimension  only).  In  general  the  quadratic  function  accounted  for 
more  of  the  variance  in  the  data  than  the  linear  one  (see  Table  4  for  a  summary  of  the 
trend  analysis  for  the  selected  pixels).  The  percentage  variance  accounted  for  ranged 
from  20.9%  for  pixel  105  to  81.9%  for  pixel  25.  The  variograms  of  selected  spectra 
were  computed  on  the  residuals  from  the  trend  -  Table  4  describes  whether  the 
residuals  from  the  linear  or  quadratic  trend  or  the  raw  data  were  used.  All  of  these 
variograms  are  periodic.  The  periodic  function  on  its  own  and  combined  with  other 
functions  such  as  power,  exponential  and  stable  exponential  functions  were 
programmed  in  GenStat  for  this  analysis.  The  stable  exponential  enables  a  concave 
shape  near  to  the  origin  to  be  fitted  satisfactorily.  The  program  is  given  in  Appendix  I. 

The  simple  periodic  function  is  given  by: 


y{h)  =  cx  cos 


2  nh 
\  CO  ) 


+  c2  sin 


r2nh 

K  (o 


where  co  is  the  wavelength,  ci  is  the  amplitude  and  C2  is  the  phase.  This  model  is  valid 
in  one  dimension  only,  which  means  that  we  can  apply  it  satisfactorily  to  the  spectra. 
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The  periodic  with  power  function  (PPF)  is  given  by: 


y{h)  -  wha  +  c,  cos; 


(  2  7ih  \ 


CO 


+  c2  sin 


2nh 


co 


\ 

y 


where  w  is  the  intensity  of  variation,  and  a  is  the  exponent  which  describes  the 
curvature. 

The  periodic  function  combined  with  an  exponential  one  (PEF)  is  given  by: 


y{h)  =  1  -  expy 


f-M 

r2nh  ^ 

( 2nh  y 

( - >  +  c,  cos 

+  c2  sin 

\r  J 

l  <°  ) 

l  ®  J 

where  r  is  the  distance  parameter  of  the  exponential  function. 

Finally  the  periodic  function  combined  with  a  stable  exponential  one  (PES)  is  given 
by: 


y(h )  =  1-expy 


-hc 


■  +  Cj  cos 


lull 


(0 


+  c2  sm 


.  f  2nh  ^ 


CO 


where  a  is  the  exponent.  This  model  takes  account  of  an  upwardly  concave  slope  in 
the  experimental  variogram  near  to  the  origin. 


This  suite  of  models  was  fitted  each  time  and  the  function  with  the  smallest  residual 
sum  of  squares  was  selected.  The  one  that  fitted  the  most  commonly  was  the  periodic 
combined  with  a  power  function.  Table  4  gives  a  summary  of  the  variogram  analysis, 
the  data  used,  the  type  of  model  fitted  and  the  model  parameters. 

Although  all  variograms  showed  evidence  of  trend  and  periodicity  -  there  are  still 
several  different  forms.  For  example  some  have  a  longer  wavelength  than  others, 
others  continue  to  increase  while  some  flatten,  some  have  larger  variances  and 
amplitudes.  The  variograms  for  pixels  1  and  25  are  similar  (as  are  their  spectra), 
Figures  12a  and  13a,  but  the  amplitudes  of  the  fitted  models  are  different,  Figures  12c 
and  13c.  This  also  applies  to  pixels  71  and  85,  Figures  15  and  22,  respectively,  where 
the  amplitudes  and  variances  are  quite  different,  which  reflect  the  differences  in  the 
spectra.  Their  variograms  have  similar  shapes,  however.  Pixels  82  and  108,  Figures 
21c  and  30c,  have  variograms  with  similar  wavelengths,  but  very  different  amplitudes 
(Table  4),  which  are  reflected  in  the  differences  in  their  spectra,  Figures  21a  and  30a. 
The  variograms  of  pixels  1  and  95,  Figures  12c  and  25c,  are  different  in  form;  the 
former  has  a  longer  wavelength,  but  both  have  similar  amplitudes  (Table  4).  Pixel 
163,  Figure  50  b  and  c,  has  a  very  different  variogram  from  all  of  the  others.  That 
computed  on  the  residuals  has  a  single  large  peak  and  is  fitted  by  a  simple  periodic 
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function  with  a  shorter  wavelength  than  those  above.  Pixel  230  has  a  distinctive 
variogram  with  a  long  wavelength  of  about  735  nm  and  a  small  amplitude,  Figure  65 
c  and  Table  4.  Pixel  250  has  a  similar  variogram  and  spectrum  to  of  pixel  82,  Figures 
68  and  21,  respectively.  Pixel  266,  Figure  73,  has  a  variogram  with  a  similar 
wavelength  to  that  of  272,  Figure  76,  but  their  amplitudes  are  different  and  there  are 
subtle  differences  in  their  spectra.  The  variograms  of  pixels  288  and  317,  Figures  81 
and  83,  respectively,  have  similar  wavelengths,  but  different  amplitudes  which  are 
reflected  by  small  differences  in  their  spectra.  Pixels  323  and  327  have  some 
differences  in  their  spectra,  although  they  are  similar  and  this  is  reflected  by  the 
different  wavelengths  of  the  variograms  (Figures  84  and  85,  and  Table  4).  The 
wavelength  for  the  former  is  776  nm  and  the  latter  254  nm,  Figures  84  and  85. 


Visual  appraisal  of  ground  cover  class  spectra  and  variograms 

The  spectra  and  variograms  for  classes  1,  4,  5  and  6  are  distinctive.  There  are  small 
differences  within  these  classes  that  are  evident  in  both  differences  in  the  spectra  and 
the  variogram  model  parameters.  The  spectra  of  class  2  are  more  heterogeneous.  The 
overall  shapes  of  the  spectra  are  similar  -  differences  relate  to  different  reflectance 
maxima  and  subtle  differences  in  the  sharpness  of  the  peaks  in  the  central  section. 
Class  5  has  similar  spectra  to  some  of  those  in  class  2  -  see  the  list  above. 


Multivariate  analysis  of  Fort  Hood  spectra 
Principal  component  analysis 

A  principal  components  analysis  was  done  using  all  of  the  pixels  and  the  spectral 
values.  Table  5  gives  the  latent  roots  and  vectors  of  the  first  three  principal 
components,  which  account  for  99.5%  of  the  variance.  Eigenvector  values  greater 
than  0.1  are  in  bold.  For  PC  1  these  are  for  wavebands  22  to  34,  and  45.  For  PC  2 
they  are  6,  48  to  49,  50  to  57,  60  to  79,  and  for  PC  3  1  to  4,  6  tol2,  49  to  57  and  82  to 
100.  Figure  10  shows  the  pixels  plotted  in  the  projection  of  PCI  and  PC  2,  using  the 
PC  scores.  There  are  four  well-defined  and  coherent  groups,  1,  3,  4  and  6  which 
accords  with  the  visual  observations  of  the  spectra  and  the  variograms.  Class  4  is  the 
most  distinct  and  this  has  spectra  with  large  average  reflectance  values  -  these  relate 
to  the  bare  ground  class  cover.  Class  2  is  diffuse  and  merges  with  class  5,  which  again 
was  apparent  from  the  visual  appraisal.  Classes  2  and  5  have  a  larger  extent  on  axis  2 
than  on  axis  1,  although  both  have  larger  extents  on  PC  1  than  any  other  class. 
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Principal  component  1 


Figure  10.  Plot  of  individual  pixels  in  the  projection  of  principal  components  1  and  2 
using  the  PC  scores.  The  ground  cover  classes  are  shown  near  to  the  relevant 
group  of  pixels,  and  the  line  shows  the  approximate  separation  between 
classes  2  and  5. 


Non-hierarchical  ft-means  dynamic  clustering 

This  is  an  alternative  approach  to  classification  compared  with  the  more  commonly 
used  hierarchical  methods.  Webster  and  Oliver  (1990)  describe  the  method.  Non- 
hierarchical  classification  is  also  known  as  dynamic  clustering.  The  population  is 
subdivided  at  a  single  level  into  as  many  classes  as  desired.  The  approach  subdivides 
a  set  of  individuals  into  two  or  more  disjoint  groups.  Each  individual  belongs  to  one, 
and  only  one ,  group.  The  general  aim  is  to  subdivide  the  population  optimally  such 
that  there  is  minimum  variation  within  the  classes  and  the  difference  between  them  is 
maximised. 
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A  mathematical  criterion  is  chosen  as  a  basis  for  optimising  the  subdivision,  such  as 
the  within  groups  sums  of  squares  or  Wilks'  criterion  (L)  to  measure  the  dispersion 
within  groups,  or  Trace  W  *B  to  measure  the  separation  between  groups.  The 
population  can  be  divided  into  an  arbitrary  number  of  groups  at  the  outset  or  the 
number  could  be  based  on  information  from  another  analysis.  The  test  criterion  is 
calculated,  and  individuals  are  then  moved  from  group  to  group  and  the  criterion 
recalculated.  If  the  change  improves  the  criterion  the  move  is  retained,  otherwise  it  is 
not.  There  are  different  ways  of  moving  individuals  from  group  to  group  in  an 
iterative  way  to  try  to  obtain  an  optimum  (Webster  and  Oliver,  1990).  In  general,  the 
optimal  number  of  groups  to  subdivide  the  population  into  is  unlikely  to  be  known, 
therefore  more  groups  than  are  likely  are  chosen  at  first.  The  number  of  groups  can 
then  be  reduced  one  at  a  time  by  fusing  the'  two  most  similar  groups  and  recalculating 
the  criterion.  Individuals  are  then  moved  as  before. 

Non-hierarchical  classification  has  an  additional  advantage  compared  with  the 
hierarchical  ones  in  that  once  individuals  are  assigned  to  groups  they  can  still  be 
moved.  This  means  that  as  the  groups  change  character  as  individuals  are  removed  or 
added  others  can  be  moved  to  optimise  the  criterion.  Once  an  individual  is  grouped  in 
a  hierarchical  method  it  is  irrevocable. 

One  of  the  difficulties  with  non-hierarchical  methods  is  that  it  is  sometimes  difficult 
to  decide  how  many  groups  is  optimal.  Groups  tend  to  be  of  a  similar  size  and  shape, 
and  weak  clusters  might  not  be  isolated  as  distinct  classes.  There  is  a  solution  to  this 
that  uses  Wilks’  criterion,  L,  and  gL  is  plotted  against  g,  where  g,  is  the  number  of 
classes. 

I  applied  this  method  clustering  to  the  pixels  from  the  Fort  Hood  image.  There  were 
too  may  wavebands  for  the  method  to  work  successfully  -  the  large  number  caused 
instability  in  the  matrices.  Therefore,  I  did  the  analysis  on  the  scores  of  the  first  three 
principal  components,  and  on  50%  and  25%  of  the  wavebands.  The  sums  of  squares 
criterion  provided  the  most  stable  results.  Starting  with  10  groups,  the  pixels  were 
reclassified  until  three  groups  remained.  It  was  not  possible  to  calculate  Wilks’ 
criterion,  therefore  I  plotted  the  sums  of  squares  criterion  against  g  for  each  data  set. 
Figure  1 1  shows  the  results  for  the  three  data  sets.  The  most  likely  number  of  groups 
is  shown  where  the  graph  goes  below  the  general  trend.  Based  on  this  the  analyses 
using  the  raw  spectra  suggest  that  five  groups  is  optimal,  whereas  for  the  PC  scores  6 
groups  is  optimal.  Figure  10  of  the  scores  plotted  in  the  plane  of  components  1  and  2 
shows  five  separate  groups.  Nevertheless,  it  is  evident  that  groups  2  and  5  although 
part  of  a  continuum,  do  occupy  distinctive  parts  of  this  projection. 

Table  6  gives  the  classes  to  which  each  pixel  belongs  for  the  three  sets  of  data.  For  the 
66  wavebands  both  the  5-group  and  6-group  classifications  are  given  for  comparison 
with  the  six  classes  from  the  PC  A  scores.  There  are  minor  differences  between  the  six 
classes  for  the  raw  data  and  the  PC  scores  in  ground  cover  class  2.  The  bold  lines 
show  the  limit  of  each  of  the  six  ground  cover  classes.  The  comparison  below  uses  the 
results  of  the  optimal  classifications  for  the  raw  data  with  66  bands  (5  groups)  and  the 
PC  scores  (6  groups). 
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Criterion 


Class  1  pixels 


For  class  1  the  classification  results  for  the  three  analyses  show  that  the  70  pixels 
belong  to  the  same  group.  The  group  is  pure  and  this  reflects  the  homogeneity  of  the 
spectra  and  the  variograms  from  a  visual  appraisal.  The  largest  reflectance  values  in 
this  class  are  about  300  nm.  The  ground  cover  for  this  class  is  grass. 

The  variograms  of  the  raw  data  are  markedly  concave  near  to  the  origin.  Figures  12b 
and  13b  of  the  raw  variograms  of  pixels  1  and  25  show  that  they  are  similar.  The 
variograms  were  computed  again  from  the  quadratic  residuals,  Figures  12c  and  13  c 
for  pixels  1  and  25,  again  they  are  very  similar.  They  were  fitted  best  by  a  periodic 
with  power  function  (PPF)  and  they  have  a  wavelength  of  about  500  nm,  and 
amplitudes  between  600  and  just  over  1000.  Table  4  gives  the  parameters  of  the  best 
fitting  models. 


Class  2  pixels 

This  is  the  largest  class  with  92  individuals  and  the  ground  cover  comprises  trees 
interspersed  with  herbaceous  ground  cover  plants.  It  does  not  appear  to  be  the  kind  of 
dense  woodland  that  occurs  at  Fort  A.  P.  Hill.  The  spectra  for  this  ground  cover  class 
are  very  different  in  shape  from  those  of  class  1,  Figures  15a  and  48a.  The 
classification  in  all  cases  shows  that  this  is  a  heterogeneous  group.  For  both  the  raw 
data  and  the  PC  scores  the  class  comprises  three  groups  from  the  non-hierarchical 
classification,  but  there  are  differences  between  them.  For  the  raw  data  the  class 
comprises  mainly  groups  3  (26  pixels)  and  4  (59  pixels),  with  only  seven  pixels 
belonging  to  group  5.  The  pixels  of  the  latter  group  are  scattered  within  the  class. 
Those  for  groups  3  and  4  are  similar  in  shape,  but  the  reflectance  values  for  group  3 
are  larger  (a  maximum  of  over  400,  compared  with  about  300  for  group  4).  The 
spectra  for  group  5  also  have  a  similar  shape  to  groups  3  and  4,  but  their  reflectance 
values  are  much  smaller,  i.e.  reaching  a  maximum  of  about  130.  The  six-group 
classification  using  PC  scores  has  6  pixels  only  in  group  3;  these  have  spectra  with  the 
largest  reflectance  values  in  this  ground  cover  class.  There  are  47  pixels  in  group  4 
(more  or  less  equivalent  to  group  3  for  the  raw  data  with  5  groups)  and  39  in  group  5 
(more  or  less  equivalent  to  group  4  for  the  raw  data  with  5  groups).  The  six-group 
classifications  for  the  raw  data  and  the  PC  scores  are  similar. 


The  experimental  variograms  of  the  raw  data  are  distinctively  different  from  those  of 
class  1,  Figures  15b  to  48b.  The  variograms  were  recomputed  on  the  residuals  and 
modelled  by  the  periodic  with  power  function.  The  model  parameters  for  selected 
pixels  are  given  in  Table  4.  All  of  the  models  have  wavelengths  of  between  465  nm 
and  476  nm.  Pixels  107  and  140  belong  to  group  3  for  the  raw  data  and  the  PC  A 
scores.  They  have  the  largest  amplitudes  of  the  variograms  in  this  class  -  over  5555, 
Figures  29  and  40.  Pixels  71,  77  and  82,  Figures  15,  17  and  21,  respectively,  belong  to 
group  3  for  the  raw  data  and  4  for  the  PC  scores.  They  have  large  amplitudes  of 
between  3669  and  4881.  Pixels  85  and  128  (Figures  22  and  36)  belong  to  group  4  for 
the  raw  data  and  the  PC  scores.  Their  amplitudes  are  intermediate  with  values  of 
about  2600.  Pixels  90,  91  and  162  (Figures  23,  24  and  48)  belong  to  group  4  for  the 
raw  data  and  group  6  for  the  PC  scores.  They  have  a  small  amplitude  of  between  1090 
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and  2197.  Finally  pixel  104  belongs  to  groupS  for  both  the  raw  data  and  PC  scores  -  it 
has  the  smallest  amplitude.  Nevertheless  the  variograms  have  a  consistent  form.  The 
differences  in  the  amplitudes  of  the  periodic  component  of  the  fitted  models  reinforce 
the  fact  that  this  group  has  differences  in  the  magnitude  of  the  reflectances. 

The  complexity  of  this  ground  cover  class  is  not  surprising  because  trees  result  in 
variable  reflectances  associated  with  the  shadows  created  by  their  crowns  and  their 
different  heights  and  extents.  In  the  1 -metre  imagery  described  in  Appendix  (II  -  the 
paper  submitted  to  the  Internationa]  Journal  of  Remote  Sensing)  these  effects  could  be 
seen  as  ‘pock-marks’. 

Class  3  pixels 

The  spectra  of  this  class  are  very  different  from  the  previous  classes.  They  are  flatter, 
and  the  reflectance  increases  gradually  from  the  small  wavelengths  to  the  longer  ones 
where  they  become  more  bulbous  in  shape.  The  spectra  have  small  maximum 
reflectance  values  of  about  130.  The  class  is  well  defined  visually  and  in  the 
classifications.  The  six-group  classifications  for  both  the  raw  data  and  the  PC  scores 
distinguish  the  class  from  pixels  in  class  2,  whereas  the  five-group  classification  of  he 
raw  data  links  these  pixels  with  the  low  reflectance  pixels  (group  5)  in  class  2.  This 
ground  cover  class  is  asphalt  and  there  should  be  no  overlap  with  class  2  spectra.  The 
group  5  pixels  in  class  2  do  not  have  a  similar  shape  to  the  spectra  of  this  class;  they 
simply  have  small  maximum  reflectance  values.  This  suggests  that  the  optimal 
classification  into  six  groups  based  on  the  PC  scores  is  more  powerful  in  terms  of 
discrimination  than  that  using  the  raw  data. 

This  group  could  also  be  discriminated  on  the  basis  that  a  linear  trend  accounted  for 
most  of  the  variance  in  the  data  (81.6%  to  89%),  whereas  for  the  majority  of  the 
spectra  most  of  the  variance  was  accounted  for  by  a  quadratic  function.  The 
experimental  variograms  of  the  raw  data  and  those  computed  from  the  linear  residuals 
are  different  from  those  for  classes  1  and  2.  The  variograms  of  the  residuals  are  a 
single  peak  and  were  fitted  best  by  a  simple  periodic  function.  Figures  49  to  64  shows 
a  selection  and  Table  4  gives  the  model  parameters.  They  have  a  wavelength  of  about 
350  nm,  and  a  small  amplitude  of  about  20. 


Class  4  pixels 

This  class  is  well  defined  both  visually  and  by  the  classifications.  It  is  pure  and  clearly 
discriminated.  The  ground  cover  is  bare  soil  and  the  maximum  reflectance  values  are 
large,  close  to  500.  The  spectra  have  an  overall  bulbous  shape  with  less  defined  peaks 
and  troughs  than  many  of  the  other  spectra. 

The  variograms  were  computed  from  the  quadratic  residuals.  Figures  65  and  66 
(pixels  230  and  234)  show  that  the  wavelengths  are  close  to  700  nm  and  the 
amplitudes  are  small,  about  150.  For  pixel  242,  Figure  67  and  Table  4,  the  wavelength 
is  smaller  -  about  450  nm. 
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Class  5  pixels 


This  ground  cover  class  is  mixed  herbaceous  vegetation.  The  spectra  have  a  similar 
general  shape  to  those  of  class  2.  For  the  five-group  classification  this  class  overlaps 
with  the  26  pixels  in  group  3,  i.e.  the  ones  with  the  larger  reflectance  values.  For  the 
six-group  classification  the  overlap  is  much  less  -  for  the  PC  scores  it  is  with  six 
pixels  only.  This  latter  seems  more  reasonable  in  that  it  is  likely  that  there  are 
herbaceous  plants  between  the  trees  in  places  in  class  2.  These  results  support  the 
statement  above  that  the  classification  using  the  PC  scores  is  more  discriminatory. 

The  variograms  of  this  class  were  computed  on  the  quadratic  residuals.  The 
wavelengths  of  the  periodic  with  power  functions  are  about  470  nm,  which  is  close  to 
those  for  class  2  for  which  the  spectra  are  also  similar.  The  raw  variograms  and  those 
on  the  residuals  have  a  similar  shapes  also,  Figures  68  to  80,  The  amplitudes  of  these 
functions  are  large,  up  to  12  OOOnm,  Table  4.  However,  there  is  a  considerable  range 
in  amplitude  from  2993  (pixel  159)  to  12  487  (pixel  284).  The  pixels  belong  to  group 
3  for  both  the  raw  data  and  PC  scores  and  have  similar  properties  to  pixels  in  this 
group  for  class  2.  This  ground  cover  class  has  the  largest  amplitudes  for  the 
variogram  models. 


Class  6  pixels 

This  groundcover  class  is  disturbed  soil,  but  it  is  a  well-defined  group  using  the  raw 
data  and  PC  scores.  Visually  these  spectra  are  also  different,  yet  the  classifications 
indicate  an  overlap  with  class  1.  The  spectra  in  these  classes  have  similar  peaks  and 
troughs,  but  they  are  less  defined  and  have  smaller  maximum  reflectance  values. 

This  ground  cover  class  has  distinctive  variograms,  as  well  as  distinctive  spectra, 
Figures  82  to  85.  The  experimental  variograms  of  the  raw  data  have  a  marked 
unbounded  form.  A  quadratic  trend  accounted  for  a  larger  percentage  of  the  variance 
than  for  any  other  class.  Variograms  were  computed  on  the  residuals  from  this  trend. 
The  best  fitting  model  was  also  different  from  the  other  classes,  being  a  periodic  with 
exponential  function  in  the  main.  The  wavelength  is  about  300  nm  on  average,  which 
is  smaller  than  those  for  the  other  classes.  The  amplitude  is  also  the  smallest  being 
about  50  on  average.  Table  4  gives  the  model  parameters  of  the  variograms.  Figures 
81  to  85  show  selected  variograms  from  the  raw  data  and  the  quadratic  residuals. 
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Summary  of  results 


For  ground  cover  classes  1,  3,  4  and  6  the  results  show  a  clear  distinction  between 
them  in  terms  of  their  spectra.  They  occupy  well-defined  parts  of  the  property  space 
in  the  projection  of  the  first  two  principal  components.  A  comparison  of  the 
classification  results  using  the  raw  data  and  the  PC  scores  suggest  that  the  PC  score 
enable  a  better  discrimination  between  the  spectra.  However,  for  the  5-group 
classification  using  the  raw  data  and  the  6-groups  classification  using  the  PC  scores 
(the  optimal  classifications)  classes  1  (grass)  and  6  (disturbed  soil),  which  have 
different  spectra,  have  been  placed  in  the  same  multivariate  group  as  each  other  by 
both  classifications.  Since  the  PC  scores  provided  the  basis  for  good  discrimination, 
they  were  classified  into  7  groups  to  see  whether  there  was  any  explanation  for  the 
overlap  between  classes  1  and  6.  At  the  7-group  level  class  6  is  now  distinct  from 
class  1.  All  of  the  other  classes  remain  similar.  It  seems  that  the  variation  in  the 
spectra  of  class  2  is  greater  than  that  between  classes  1  and  6.  Therefore,  caution  is 
required  in  taking  the  prescribed  route  through  the  analysis.  The  other  puzzle  is  the 
overlap  between  classes  2  and  5  which  possibly  relates  to  the  fact  that  there  is  ground 
vegetation  between  the  trees. 

In  spite  of  some  anomalies  the  analyses  suggest  that  a  combination  of  visual 
inspection,  multivariate  classification  and  ordination,  and  variography  enable  us  to 
distinguish  between  different  kinds  of  ground  cover  class  in  terms  of  their 
hyperspectral  information.  The  spectra  of  classes  1,  3,  4  and  6  were  clearly  different. 
The  complexity  of  class  2  spectra  can  probably  be  explained  by  the  texture  in  the 
pixels  arising  from  different  tree  heights,  extents  and  proximity,  with  areas  of  ground 
cover  vegetation  interspersed  and  showing  through  in  parts. 
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b) 


Figure  17.  Ft  Hood:  pixel  77:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals. 
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Figure  18.  Ft  Hood:  pixel  78:  a)  spectrum,  b)  experimental  variogram  of 
quadratic  residuals. 
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Figure  21.  Ft  Hood:  pixel  82:  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals, 
a) 


b) 


Figure  22.  Ft  Hood:  pixel  85:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals. 
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Figure  23.  Ft  Hood:  pixel  90:  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals. 
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Figure  24.  Ft  Hood:  pixel  91:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals. 
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Figure  25.  Ft  Hood:  pixel  95:  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals. 
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Figure  26.  Ft  Hood:  pixel  96:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals. 
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a) 


and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals. 
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Figure  28.  Ft  Hood:  pixel  104:  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals. 
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a) 


Figure  29.  Ft  Hood:  pixel  107:  a)  spectrum,  b)  variogram  and  fitted  model  (PEF)  of 
raw  values  clipped,  c)  variogram  and  fitted  model  (PPF)  of  quadratic  residuals, 
a) 


b) 


Figure  30.  Ft  Hood:  pixel  108:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals. 
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Figure  3 1  Ft  Hood:  pixel  111:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals, 
a) 


b) 


Figure  32  Ft  Hood:  pixel  114:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals. 
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and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals, 
a) 


b) 


Figure  34.  Ft  Hood:  pixel  121:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals. 
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Figure  36.  Ft  Hood:  pixel  128  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals, 
a) 


b) 


Figure  37.  Ft  Hood:  pixel  131  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals. 
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a) 


and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals, 
a) 


and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals. 
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b) 


Figure  40.  Ft  Hood:  pixel  140:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals, 
a) 


b) 


Figure  41 .  Ft  Hood:  pixel  141  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals. 
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Figure  42.  Ft  Hood:  pixel  142:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals, 
a) 
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Figure  43.  Ft  Hood:  pixel  143:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals. 
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Figure  44.  Ft  Hood:  pixel  140:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals, 
a) 
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c) 


Figure  45.  Ft  Hood:  pixel  159:  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals. 
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Figure  48.  Ft  Hood:  pixel  162:  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals. 
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Figure  49.  Ft  Hood:  pixel  163:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  raw  values,  and  c)  as  b)  but  clipped  at  500  nm  lag  distance, 
a) 
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Figure  50.  Ft  Hood:  pixel  163  a)  spectrum,  b)  variogram  and  fitted  model  (PPF)  of 
raw  values,  c)  variogram  and  model  (periodic)  of  linear  residuals,  clipped. 
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Figure  51  Ft  Hood:  pixel  167:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  raw  values, 
a) 
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Figure  52.  Ft  Hood:  pixel  168:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  raw  values. 
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a) 


Figure  53.  Ft  Hood:  pixel  169:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  raw  values, 
a) 


b) 


Figure  54.  Ft  Hood:  pixel  171:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  raw  values. 
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Figure  55.  Ft  Hood:  pixel  188:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  raw  values, 
a) 


b) 


Figure  56.  Ft  Hood:  pixel  190:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  raw  values. 


a) 


b) 


Figure  61 .  Ft  Hood:  pixel  206:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  raw  values. 


a) 


Wavelength  nm 


b) 


Figure  62.  Ft  Hood:  pixel  207:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  raw  values,  c)  as  b)  but  clipped  at500  nm  lag  distance. 
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a) 


b) 


Figure  64.  Ft  Hood:  pixel  214:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  raw  values. 
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a) 


model  (PES)  of  raw  values,  c)  experimental  variogram  and  fitted  model  (PEF) 
of  raw  values,  clipped  at  500nm  lag  distance,  d)  experimental  variogram  and 
fitted  model  (PEF)  of  quadratic  residuals  and  e)  experimental  variogram  and 
fitted  model  (PPF)  of  quadratic  residuals,  clipped  at  300  nm  lag  distance. 
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Figure  68.  Ft  Hood:  pixel  250:  a)  spectrum,  b)  experimental  variogram  of  raw  values 
and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals, 
a) 
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Figure  69.  Ft  Hood:  pixel  255:  a)  spectrum,  b)  experimental  variogram  of  raw  values 
and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals. 


b) 


Figure  70.  Ft  Hood:  pixel  257:  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals. 

a) 


b) 


Figure  71.  Ft  Hood:  pixel  260:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PPF)  of  quadratic  residuals. 
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Figure  74. 
a) 


Ft  Hood:  pixel  268:  a)  spectrum,  b)  experimental  variogram  of 
quadratic  residuals. 
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Figure  75.  Ft  Hood:  pixel  269:  a)  spectrum,  b)  experimental  variogram  of 
quadratic  residuals. 


Lag  distance  nm 

Figure  76.  Ft  Hood:  pixel  272:  a)  spectrum,  b)  experimental  variogram  of  raw  values 
and  c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals, 
a) 


b) 


Figure  77.  Ft  Hood:  pixel  273:  a)  spectrum,  b)  experimental  variogram  of 
quadratic  residuals. 


Wavelength  nm 


Figure  78. 
a) 


Ft  Hood:  pixel  278:  a)  spectrum,  b)  experimental  variogram  of  raw 
values  and  c)  experimental  variogram  of  quadratic  residuals. 


Figure  79.  Ft  Hood:  pixel  282:  a)  spectrum,  b)  experimental  variogram  of  raw 
values  and  c)  experimental  variogram  of  quadratic  residuals. 
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Figure  81.  Ft  Hood:  pixel  288:  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
c)  experimental  variogram  and  fitted  model  (PEF)  of  quadratic  residuals,  d) 
variogram  and  fitted  model  (PEF)  of  quadratic  residuals  clipped  at  500  nm. 


a) 


d) 


Figure  82.  Ft  Hood:  pixel  300:  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
c)  experimental  variogram  and  fitted  model  (PEF)  of  quadratic  residuals,  d) 
variogram  and  fitted  model  (PES)  of  quadratic  residuals  clipped  at  500  nm,  e) 
variogram  and  fitted  long  wavelength  model  (PPF)  of  quadratic  residuals. 
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Figure  83.  Ft  Hood:  pixel  317:  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
c)  experimental  variogram  and  fitted  model  (PPF)  of  quadratic  residuals,  d) 
variogram  and  fitted  model  (PEF)  of  quadratic  residuals  clipped  at  500  nm. 
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Wavelength  nm 
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Figure  84.  Ft  Hood:  pixel  323:  a)  spectrum,  b)  experimental  variogram  and  fitted 
long  wavelength  model  (PPF)  of  quadratic  residuals. 

a) 


Figure  85.  Ft  Hood:  pixel  327:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  (PEF)  of  quadratic  residuals,  c)  as  b)  fitted  to  long  wavelength  model  (PEF). 


Table  4:  Summary  of  trend  analysis  and  model  parameters  of  variograms  for  selected  pixels  for 
Fort  Hood 


%  trend  removed 

linear  quadratic 

6. 


Trend  Model  Fitted 

Used  PPF  PES  PEF  Nugget 


80.6  Q 


81.9  Q 


53.0  Q 


55.5  Q 


53.0  Q 


51.4  Q 


51.6  Q 


51.4  Q 


57.6  Q 


52.4  Q 


51.6  Q 


51.5  Q 


52.3  Q 


Q  F 

None  F 


54.1  Q 


52.4  Q 


53.9  Q 


67.1  Q 


51.1  Q 


53.5  Q 


55.4  Q 


62.4  Q 


58.0  Q 


49.8  Q 


54.4  Q 


53.9  Q 


53.6  Q 


59.0  Q 


53.2  Q 


88 

88 


89.3 

89.3 


87.0  None 


86.9  None 


Model  Parameters 

get  Sill 

Distance  Gradient  Exponent 

Wavelength 

Amplitude  Phase 

129.3  0.4205 

527.6 

605.6 

2.747 

1808 

130.6 

473.7 

525.1 

3.090 

1194  0.4399 

472.4 

4881 

2.358 

F 


F 


F 

F 


Periodic 

F 


886.7  0.4514 


821.5  0.4743 


497.4  0.5041 


428.8  0.5000 


193.5  0.5263 


90.74  0.5411 


375.3  0.4954 


190.8  0.5240 


48.82  0.5615 


1425  0.4626 


108 

3.7 

56.3 

Q 

F 

0 

82.99 

0.5672 

467.8 

509.4 

2.332 

111 

1.1 

53.1 

Q 

F 

0 

354.2 

0.4997 

471.0 

1794 

2.301 

114 

5.2 

54.7 

Q 

F 

58.79 

0.5397 

466.0 

331.5 

2.327 

120 

3.0 

56.5 

Q 

F 

84.09 

0.5860 

463.7 

569.9 

2.269 

217.7  0.5352 


87.69  0.5184 


519.1  0.4939 


356.4  0.5052 


535.2  0.4735 


1032  0.4645 


1  1241 

0.4691 

476.2 

5555 

2.381 

1  271.5 

0.5555 

465.9 

1621 

2.345 

0 

271.3 

0.5439 

467.0 

1587 

2.305 

0 

73.55 

0.5527 

465.5 

460.1 

2.218 

0 

65.16 

0.5997 

465.8 

461.6 

2.228 

0 

902.2 

0.4607 

467.0 

3934 

2.333 

0 

290.2 

0.4886 

469.8 

1346 

2.401 

0 

40.68 

0.5993 

447.6 

281.5 

2.173 

0 

293.7 

0.5126 

465.1 

1570 

2.259 

0 

0.1532 

1.0930 

385.2 

20.38 

2.199 

0 

0.1178 

1.1400 

400.8 

24.34 

2.12 

0 

39.15 

2.654 

0 

0.098 

1.1820 

22.86 

2.266 

0 

0.1983 

1.0660 

399.3 

22.42 

2.141 

0 

0.462 

0.9367 

422.5 

23.77 

2.136 

Q  quadratic  detrend 
Lin  linear  detrend 


F  is  best  fitting  model 
PPF  periodic  with  power  function 
PES  periodic  with  stable  function 
PEF  periodic  exponential  function 


L  fitted  model  with  longer  wavelength 
C  fitted  model  with  shorter  lag  distance 


Table  4:  Summary  of  trend  analysis  and  model  parameters  of  variograms  for  selected  pixels  for 
Fort  Hood  (continued) 


Pixel 

%  trend  removed  i 

linear  Quadratic! 

Trend 

Used 

Model  Fitted 

PPF  PES  PEF 

Nugget 

Sill 

Model  Parameters 

Distance  Gradient  Exponent  Wavelength  Amplitude  Phase 

169 

86.1 

86.9 

None 

F 

0 

0.2001 

1 .0480 

400.3 

23.16 

3.145 

171 

83.1 

83.3 

None 

F 

0 

0.2163 

1.0080 

390.3 

21.54 

2.264 

188 

82.0 

81.9 

None 

F 

0 

0.0923 

1.1070 

385.4 

18.78 

2.171 

190 

81.6 

82.1 

None 

F 

0 

0.5801 

0.8736 

400.8 

22.47 

2.254 

193 

84.8 

84.8 

None 

F 

0 

0.0311 

1.2780 

385.6 

18.15 

2.113 

195 

85.5 

86.8 

None 

F 

0 

0.4117 

0.9324 

404.8 

20.78 

2.167 

198 

82.3 

83.8 

None 

F 

0 

0.8126 

0.7941 

412.2 

21.71 

2.256 

204 

88.0 

89.2 

None 

F 

0 

0.1689 

1.1010 

381.9 

19.52 

2.192 

206 

87.6 

None 

F 

0 

0.0761 

1.1880 

370.7 

19.81 

2.226 

207 

84.1 

84.2 

None 

F 

0 

0.1049 

1.0910 

19.93 

2.181 

II 

84.1 

84.2 

None 

B 

0 

0.0193 

1.3890 

429.3 

26.52 

1.985 

213 

82.5 

83.0 

None 

mmm 

0 

0.2171 

0.9736 

410.8 

18.96 

2.166 

214 

87.9 

88.8 

None 

F 

0 

0.1531 

1.1420 

367.1 

21.61 

2.315 

230 

60.4 

97.0 

Q 

F 

0 

409.8 

52.89 

734.6 

150.7 

1.855 

234 

61.5 

96.7 

Q 

F 

0 

449.7 

58.03 

669.5 

138.4 

2.063 

ll 

61.5 

96.7 

None 

F 

0 

12719 

348.30 

373.9 

1775 

1.637 

242 

60.8 

97.7 

Q 

F 

0 

336.6 

39.34 

411.5 

65.93 

1.643 

'• 

60.8 

97.7 

Q 

F 

0 

37.07 

0.4111 

106.1 

54.13 

1.892 

ll 

60.8 

97.7 

None 

F 

0 

15395 

423.00 

1.1020 

367.2 

1733 

1.169 

li 

60.8 

97.7 

None 

F 

0 

7804 

126.90 

476.3 

2860 

1.649 

250 

54.5 

Q 

F 

0 

1228 

0.4397 

472.7 

4823 

2.409 

255 

■NiB 

51.7 

Q 

F 

0 

2919 

0.3735 

483.5 

9268 

2.464 

257 

■BB 

52.6 

Q 

F 

0 

2814 

0.4079 

477.8 

10054 

2.433 

260 

0.2 

54.6 

Q 

F 

0 

956.1 

0.4575 

471.2 

3993 

2.383 

261 

0.1 

55.3 

Q 

F 

0 

578.8 

0.4572 

470.0 

2420 

2.398 

ifl 

0.0 

54.1 

Q 

F 

0 

1984 

0.4187 

476.3 

7277 

2.436 

trn 

0.2 

50.2 

Q 

F 

0 

2885 

0.4289 

475.1 

11555 

2.326 

274 

0.0 

57.6 

Q 

F 

0 

1315 

0.4071 

480.5 

4497 

2.555 

284 

1.2 

48.4 

Q 

F 

0 

2848 

0.4475 

472.8 

12487 

2.268 

288 

93.0 

Q 

F 

0 

311.3 

83.19 

250.1 

65.78 

1.640 

II 

93.0 

Q 

F 

0 

325.0 

100.30 

349.2 

52.31 

2.712 

II 

93.0 

Q 

F 

0 

0.132 

1.2890 

759.6 

198.5 

1.647 

300 

40.8 

95.9 

Q 

F 

0 

210.3 

78.43 

248.5 

46.78 

1.436 

ll 

40.8 

95.9 

Q 

F 

0 

213.0 

81.98 

1.578 

274.4 

37.62 

1.373 

ll 

40.8 

95.9 

Q 

F 

0 

0.9401 

0.9040 

724.8 

108.8 

1.749 

317 

52.3 

96.9 

Q 

F 

0 

148.1 

53.86 

269.1 

37 

1.748 

ll 

52.3 

96.9 

Q 

F 

0 

24.68 

0.3342 

429.2 

60.6 

3.122 

ll 

52.3 

96.9 

Q 

F 

0 

164.4 

139.40 

697.1 

64.28 

1.962 

323 

57.2 

84.7 

Q 

F 

0 

0.037 

1.6480 

776.1 

491.8 

1.757 

C 


L 


C 

C 


c 

L 

c 

L 

C 

L 


Q  quadratic  detrend 
Lin  linear  detrend 


F  is  best  fitting  model 
PPF  periodic  with  power  function 

PES  periodic  with  stable  function  L  fitted  model  with  longer  wavelength 

PEF  periodic  exponential  function  C  fitted  model  with  shorter  lag  distance 


a 


Table  4:  Summary  of  trend  analysis  and  model  parameters  of  variograms  for  selected  pixels  for 
Fort  Hood  (continued) 


Pixel 

%  trend  removed 

Trend 

Model  Fitted 

Model  Parameters 

linear  quadratic 

Used 

PPF  PES  PEF 

Nugget 

Sill 

Distance  Gradient  Exponent  Wavelength  Amplitude  Phase 

326 

55.5 

97.2 

Q 

F 

0 

147.5 

59.83 

271.6 

31.56 

1.782 

« 

55.5 

97.2 

Q 

F 

0 

199.0 

104.60 

272.2 

37.52 

1.604 

II 

55.5 

97.2 

Q 

F 

0 

164.8 

158.80 

758.3 

65.12 

1.900 

327 

52.6 

97.3 

Q 

F 

0 

125.5 

52.21 

253.5 

31 .74 

1.755 

II 

52.6 

97.3 

Q 

F 

0 

140.3 

130.10 

673.4 

47.6 

1.972 

329 

F 

282.7 

75.50 

250.8 

65.93 

1.650 

Q  quadratic  detrend 
Lin  linear  detrend 


F  is  best  fitting  model 
PPF  periodic  with  power  function 

PES  periodic  with  stable  function  L  fitted  model  with  longer  wavelength 

PEF  periodic  exponential  function  C  fitted  model  with  shorter  lag  distance 


Table  5.  Eigenvalues  and  eigenvectors  for  Fort  Hood  using  all  spectra 
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Table  5.  Continued. 
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Table  6.  Classes  from  the  non-hierarchical  classification  of  Fort  Hood  spectra  using  66  wavebands  and  the  first  thee  principal  components 
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Fort  A.  P.  Hill 


Pixels  from  a  hyperspectral  image  of  Fort  A.  P.  Hill  (Virginia)  were  selected 
according  to  different  types  of  ground  cover  comprising  different  types  of  vegetation, 
trees,  soil  and  water.  The  image  is  shown  in  Figure  86a.  The  six  groups  of  pixels  were 
contiguous  and  relate  to  the  areas  marked  by  the  arrows.  Figure  86b  shows  the 
relative  sizes  of  the  classes  of  pixels  and  their  locations. 


Deciduous  woodland 


Less  dense  woodland 


Pine  plantation 
woodland 


Water 


Soil  with  vegetation 


Grass 


Figure  86.  a)  Colour  image  of  part  of  Fort  A.  P.  Hill,  and  b)  the  locations  and  sizes 
of  the  classes  of  pixels. 

There  were  126  wavelength  values  (nm)  in  each  spectrum  and  a  total  of  172  pixels  of 
about  3  m  x  3  m.  Table  7  gives  the  number  of  sites  in  each  class  and  the  associated 
ground  cover  class. 


Table  7.  Classes  of  different  ground  cover  type  and  number  of  pixels  in  each  for  Fort 
A.  P.  Hill. 


Class 

Type  of  ground 
cover 

Number  of 

pixels 

Pixel  range 

1 

Soil  with 

vegetation 

43 

1-  43 

2 

Plantation 

woodland? 

36 

44-  79 

3 

Shrubs? 

24 

80-103 

4 

Mixed  woodland? 

35 

104-138 

5 

Grass 

21 

139-159 

6 

Water 

13 

160-172 
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As  for  Fort  Hood  the  spectrum  for  each  pixel  was  plotted  as  a  trace  and  examined  for 
differences  and  similarities.  Experimental  variograms  were  computed  from  the 
spectral  data  for  each  pixel.  The  full  set  of  spectra  and  variograms  can  be  provided, 
but  to  simplify  interpretation  we  focus  on  selected  ones  with  particular  properties. 
This  selection  was  based  on  the  form  of  the  spectrum  and  its  variogram.  This  part  of 
the  investigation  was  done  without  reference  to  the  classes  of  ground  cover  in  case 
there  were  differences  in  the  spectra  for  a  given  classes  which  might  be  overlooked  by 
focusing  on  the  classes  only. 


Spectra 

Most  of  the  spectra  for  Fort  A.  P.  Hill  have  a  similar  form  -  the  main  differences 
relate  to  the  fine  detail  and  in  particular  their  maximum  reflectance  values.  The 
spectra  for  class  6  (water)  are  distinctly  different  from  the  others  and  from  any  of  the 
spectra  of  Fort  Hood.  Selected  spectra  are  shown  in  Figures  89  to  1 15. 

A  visual  inspection  of  the  spectra  suggests  that  1  to  159  have  a  similar  overall  shape. 
Some,  such  as  pixel  40,  Figure  91,  have  a  larger  spike  at  about  1300  nm  than  others 
such  as  pixel  4,  Figure  90.  The  other  main  difference  apart  from  small  details  in  the 
shape  is  the  average  reflectance  value.  For  class  1  the  average  spectral  value  is  about 
1800,  for  class  2  it  is  about  1 100,  for  class  3  it  isl400,  for  class  4  which  has  the  widest 
range  in  mean  values  for  the  spectrum  it  is  about  1700,  for  class  5  it  is  about  2100, 
and  for  class  6  it  is  160.  Thus  there  is  a  considerable  range  in  reflectance  values  and 
apart  from  class  6  it  is  likely  to  be  the  main  basis  for  discrimination. 


Variography 

All  of  the  experimental  variograms  showed  weak  evidence  of  trend  with  their 
unbounded  form.  Linear  and  quadratic  functions  were  fitted  to  the  wavelengths  (one 
dimension  only).  The  quadratic  function  accounted  for  more  of  the  variance  in  the 
data  than  the  linear  one  (see  Table  8  for  a  summary  for  the  selected  pixels).  It  ranged 
from  41%  for  pixel  107  to  59%  for  pixel  160.  The  variograms  of  selected  spectra  were 
computed  on  the  residuals  from  quadratic  trend  function.  For  pixels  1-159  they  have  a 
similar  form.  It  was  possible  to  fit  a  periodic  function  with  a  short  wavelength  as  well 
as  one  with  a  longer  wavelength  to  all  of  the  variograms,  Figure  89  for  example.  Most 
of  the  models  with  a  short  wavelength  were  the  periodic  with  exponential  or  stable 
exponential  function  (PEF  and  PES,  respectively).  In  general  the  shorter  wavelength 
is  about  250  nm,  and  the  longer  one  varies  between  520  nm  to  790  nm.  The  latter 
were  generally  associated  with  the  periodic  with  power  function  (PPF).  The  main 
difference  between  the  variograms  of  the  different  classes  of  ground  cover  relates  to 
the  amplitudes  of  the  fitted  models. 
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Multivariate  analysis  of  Fort  A.  P.  Hill  spectra 
Principal  component  analysis 

A  principal  components  analysis  was  done  using  all  of  the  pixels  and  the  spectral 
values.  Table  9  gives  the  latent  roots  and  vectors  of  the  first  three  principal 
components,  which  account  for  96.4%  of  the  variance.  Large  eigenvector  values  are  in 
bold.  For  PC  1  (•  0.0975)  these  are  for  wavebands  70  to  93.  For  PC  2  (•  0.13)  they 
are  for  11  to  17,  and  22  to  33,  and  for  PC  3  (•  0.4)  wavebands  63  to  65.  Figure  87 
shows  the  pixels  plotted  in  the  projection  of  PCI  and  PC  2  using  the  PC  scores.  There 
are  three  well-defined  and  coherent  groups,  1,  5  and  6.  Class  6  is  the  most  distinct  and 
this  has  spectra  with  the  smallest  average  reflectance  values.  Classes  2,  3  and  4  merge 
in  the  centre  of  the  figure.  Nevertheless,  based  on  the  PC  scores,  they  occupy  a 
distinct  part  of  the  space. 


Figure  87.  Plot  of  individual  pixels  in  the  projection  of  principal  components  1  and  2 
using  the  PC  scores.  The  ground  cover  classes  are  shown  near  to  the  relevant 
group  of  pixels,  and  the  lines  shows  the  approximate  separation  between 
classes  2,  3  and  4. 
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Non-hierarchical  fc-means  dynamic  clustering 

I  applied  this  method  clustering  to  the  pixels  from  the  Fort  A.  P.  Hill  image.  There 
were  too  may  wavebands  for  the  method  to  work  successfully  as  before.  Therefore,  I 
did  the  analysis  on  the  scores  of  the  first  three  principal  components,  and  on  63 
(alternate)  wavebands.  The  sums  of  squares  criterion  provided  the  most  stable  results. 
Starting  with  10  groups,  the  pixels  were  reclassified  until  three  groups  remained.  As  it 
was  not  possible  to  calculate  Wilks’  criterion,  I  plotted  the  sums  of  squares  criterion 
against  g  for  each  data  set.  Figure  88  shows  the  results  for  the  two  data  sets.  The  most 
likely  number  of  groups  is  shown  where  the  graph  goes  below  the  general  trend. 
Based  on  this  the  analysis  using  the  raw  spectra  the  graph  suggests  that  four  or  seven 
groups  is  optimal,  whereas  for  the  PC  scores  6  groups  is  optimal.  Figure  87  of  the 
scores  plotted  in  the  plane  of  components  1  and  2  shows  four  separate  groups. 
Nevertheless,  it  is  evident  that  groups  2  to  4  form  a  continuum  in  this  projection. 


a)  b) 


Figure  88  Sum  of  squares  criterion  x  g2  plotted  against  g  (number  of  groups)  for: 

a)  63  wavebands,  and  b)  the  scores  of  the  first  five  principal  components 


Table  10  gives  the  classes  to  which  each  pixel  belongs  for  the  two  sets  of  data.  For  the 
63  wavebands  both  the  6-group  and  7-group  classifications  are  given  for  comparison 
with  the  six  classes  from  the  PC  A  scores.  In  general  the  7 -group  classification  does 
not  improve  the  discrimination  between  classes,  therefore  it  is  not  discussed  below  in 
detail.  The  bold  lines  show  the  extent  of  each  of  the  six  ground  cover  classes. 
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Class  1  pixels 


Ground  cover  class  1  is  soil  with  vegetation.  For  the  PCA  6-group  classification  based 
on  the  scores  it  is  a  pure  group.  For  the  raw  spectra  six  pixels  belong  to  another  group 
that  overlaps  with  class  4.  This  class  occupies  a  distinct  part  of  the  property  space 
projection,  Figure  85  It  has  the  second  largest  reflectance  values.  Figures  89  to  91 
show  the  spectra  for  pixels  1 , 4  and  40. 

The  variograms  of  pixels  1  and  40  have  a  wavelength  of  about  250  nm  for  the  shorter 
wavelength  and  over  700  nm  for  the  longer  periodic  function.  The  amplitudes  of  the 
models  for  this  class  are  the  second  largest,  Table  8. 

Class  2  pixels 

This  ground  cover  class  appears  to  be  dense  pine  plantation  woodland.  For  the 
classification  based  on  PC  scores  six  pixels  belong  to  group  6  and  the  rest  to  group  4 
of  the  multivariate  classification.  For  the  raw  spectra  the  class  is  less  pure  having 
three  groups  present  (six  pixels  belong  to  group  4,  five  to  group  6  and  the  rest  to 
group  5).  Figures  92  to  99  show  the  spectra  selected  to  represent  this  class.  There  is 
more  variation  present  than  in  class  1 . 

The  wavelengths  of  the  variogram  models  fitted  with  the  longer  period  range  between 
about  558  nm  and  750  nm,  Table  8.  The  amplitudes  are  intermediate  in  terms  of  the 
full  data,  but  it  is  evident  that  as  the  wavelength  increases  so  does  the  amplitude, 
Table  8. 

Class  3  pixels 

This  is  the  least  pure  class  for  the  classification  based  on  PC  scores.  There  are  three 
multivariate  groups  in  this  class:  six  pixels  belong  to  group  4,  one  to  group  5  and  the 
rest  to  group  6.  For  the  raw  spectra  the  distribution  of  classes  is  similar  to  that  for  the 
PC  scores.  In  property  space  this  ground  cover  class  is  transitional  to  classes  2  and  4. 
It  appears  to  be  less  dense  vegetation  than  that  of  class  2  -  it  is  possibly  young  pine 
plantation.  Figures  100  and  101  are  examples  of  the  spectra  of  this  class. 

The  variograms  have  been  fitted  by  models  with  a  short  and  a  long  wavelength, 
Figures  100  and  101.  The  former  is  about  250  nm  and  the  longer  one  is  750  nm.  These 
variograms  have  a  larger  amplitude  than  those  of  class  2,  Table  8. 


Class  4  pixels 

This  ground  cover  class  is  deciduous  woodland.  For  the  PC  scores  the  class  is  almost 
pure  with  just  two  pixels  belonging  to  another  group.  For  the  raw  spectra  there  are 
three  groups  present,  with  13  pixels  belonging  to  groups  other  than  the  dominant  one. 
Figures  102  to  108  are  examples  of  spectra  belonging  to  this  group.  Their  overall 
shapes  are  similar,  but  the  height  of  the  peak  at  1400  nm  varies. 

The  variograms  of  this  class  have  wavelengths  for  the  longer  periodic  function  of 
between  522  nm  and  789  nm.  Their  amplitudes  are  larger  than  for  classes  2  and  3, 
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Table  8.  The  approximate  wavelength  of  the  shorter  fitting  periodic  function  is  about 
250  nm. 


Class  5  pixels 

This  ground  cover  class  is  grass  and  it  is  clearly  distinct  in  the  property  space 
projection  from  all  of  the  other  classes.  For  the  PC  scores  the  group  is  pure.  It  is 
reasonably  coherent  for  the  raw  spectra,  however,  the  majority  of  pixels  belong  to  the 
same  group  as  the  majority  of  pixels  in  class  1.  There  is  some  similarity  between  these 
ground  cover  classes  in  the  image,  Figure  1.  However,  the  intensity  of  vegetation 
cover  is  less  in  class  1 . 

Figures  109  and  110  are  examples  of  the  spectra  for  this  class.  The  variograms  in  this 
class  Figures  109  and  1 10  and  Table  8  have  the  largest  amplitudes  and  wavelengths  of 
between  615  nm  and  714  nm.  There  is  also  a  shorter  fitting  periodic  function  with  a 
wavelength  of  about  250  nm. 


Class  6  pixels 

This  is  the  most  distinctive  class  in  terms  of  the  raw  spectra,  Figureslll  to  168  ,  and 
the  variograms.  It  is  water  and  the  reflectance  values  are  the  smallest.  The 
multivariate  classifications  have  resulted  in  pure  classes  for  both  the  PC  scores  and 
the  raw  spectra. 

The  variograms  have  a  large  wavelength  (for  the  longer  fitting  periodic  function)  and 
the  smallest  amplitude,  Figures  111  to  168  and  Table  8.  The  wavelength  of  the  shorter 
fitting  function  is  about  250  nm. 


Summary  of  results 

The  classification  based  on  the  PC  scores  has  produced  more  coherent  and  pure 
classes  than  that  using  the  raw  data.  The  spectra  for  Fort  A.  P.  Hill  are  more  similar  in 
form  for  classes  1  to  5  than  those  for  Fort  Hood.  Nevertheless  the  classes  are 
distinguishable  by  multivariate  classification  using  the  PC  scores.  The  data  reduction 
to  three  components  compared  with  half  of  the  spectra  clearly  refines  the  analysis. 
The  shapes  of  the  variograms  for  these  data  are  less  distinctive  between  ground  cover 
classes,  but  the  amplitudes  of  the  periodic  functions  could  be  used  for  discrimination. 
The  variogram  for  the  water  is  distinguishable  from  all  of  the  others. 
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Conclusions 


This  project  has  analysed  pixel  information  for  selected  classes  of  ground  cover  in 
two  ways.  Part  II  of  the  report  analysed  a  set  of  pixels  belonging  to  a  coherent  ground 
cover  class  using  NIR.  The  results  show  that  the  classes  had  distinguishable 
characteristics  that  could  be  identified  from  the  structure  of  the  variogram.  The 
suggestion  is  that  the  texture  of  a  ground  cover  class  in  an  image  has  an  identifiable 
form. 

Part  III  of  the  report  describes  the  investigation  of  the  hyperspectral  data  for 
individual  pixels.  The  results  show  that  certain  ground  cover  classes  have  very 
distinctive  spectra,  for  instance  at  Fort  Hood.  However  when  the  ground  cover  classes 
are  different  kinds  of  trees  the  distinction'  is  less  clear  from  the  spectra  alone.  The 
variograms  provide  clues  as  to  the  differences  in  terms  of  the  sizes  of  the  variances,  in 
particular.  However  for  some  classes  the  spectra  have  distinct  wavelengths  for  the 
periodic  functions.  The  principal  component  analyses  have  been  particularly  effective 
in  concentrating  the  vital  information  and  in  so  doing  enable  better  discrimination 
between  the  ground  cover  classes  in  terms  of  the  spectral  values.  The  non-hierarchical 
multivariate  classifications  using  the  PC  scores  are  powerful  in  distinguishing  the 
classes. 

The  results  suggest  that  PCA  with  the  non-hierarchical  classification  provides  a  rapid 
means  of  discrimination.  The  variograms,  although  distinctive  for  certain  ground 
cover  types,  do  not  provide  a  rapid  means  of  discrimination  because  the  differences 
are  often  more  subtle.  For  example  at  Fort  A.  P.  Hill  for  classes  1  to  5  the  difference 
in  the  variograms  according  to  the  classes  was  the  amplitudes  of  the  fitted  periodic 
functions. 
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Lag  distance  nm 


Lag  distance  nm 


Figure  89.  A.P.  Hill:  pixel  1:  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
c)  and  d)  experimental  variogram  and  fitted  models  of  quadratic 
residuals  -  short  wavelength  PES  and  long  wavelength  PPF, 
respectively. 
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Figure  90.  A.P.  Hill:  pixel  4:  a)  spectrum,  b)  experimental  variogram  of  raw  values 
and  c)  experimental  variogram  of  quadratic  residuals. 


b) 


fitted  models  of  quadratic  residuals  -  short  and  long  wavelength  PPF  model. 
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Figure  92.  A.P.  Hill:  pixel  44:  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
c)  experimental  variogram  of  linear  residuals,  d)  and  e)  experimental 
variogram  and  fitted  models  of  quadratic  residuals  -  short  wavelength 
PEF  and  long  wavelength  PPF,  respectively. 
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a) 


Figure  93.  A.P.  Hill:  pixel  50:  a)  spectrum,  b)  &  c)  experimental  variogram  and  fitted 
models  of  quadratic  residuals  -  short  PES  and  long  wavelength  PPF  model. 

a) 


Figure  94.  A.P.  Hill:  pixel  54:  a)  spectrum,  b)  &  c)  experimental  variogram  and  fitted 
models  of  quadratic  residuals  -  short  (PES)  and  long  wavelength  (PPF  model). 
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a) 


fitted  models  of  quadratic  residuals  -  short  PES  and  long  wavelength  PPF 
model. 


a) 
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b)  c)  d) 

Figure  96.  A.P.  Hill:  pixel  70:  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
c)  variogram  of  linear  residuals  and  d)  variogram  of  quadratic  residuals. 
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a) 


b) 


Figure  97.  A.P.  Hill:  pixel  75:  a)  spectrum,  b)  experimental  variogram  and  fitted 
model  of  quadratic  residuals  (short  wavelength  PES  model). 

a) 
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Figure  98.  A.P.  Hill:  pixel  76:  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
c)  and  d)  variograms  and  fitted  models  of  quadratic  residuals  -  short 
wavelength  PES  model  and  long  wavelength  PPF  model. 
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Wavelength  nm 


Figure  99.  A.P.  Hill:  pixel  79:  a)  spectrum,  b)  experimental  variogram  of  raw  values 
c)  and  d)  variograms  and  fitted  models  of  quadratic  residuals  -  short 
wavelength  PES  model  and  long  wavelength  PPF  model. 


Figure  100.  A.P.  Hill:  pixel  87  a)  spectrum,  b)  experimental  variogram  of  raw  values, 
c)  experimental  variogram  of  linear  residuals,  d)  and  e)  experimental 
variogram  and  fitted  models  of  quadratic  residuals  -  short  wavelength 
PES  and  long  wavelength  PPF,  respectively. 
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a) 
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Figure  103.  A.P.  Hill:  pixel  107:  a)  spectrum,  b)  variogram  of  raw  values,  c)  &  d) 
variogram  and  fitted  models  of  quadratic  residuals  -  short  PES  and  long  wave  PPF. 
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a) 


b) 


Figure  104.  A.P.  Hill:  pixel  108:  a)  spectrum,  b)  experimental  variogram  of  quadratic 
residuals. 


a) 


Figure  105  .A.  P.  Hill:  pixel  115:  a)  spectrum,  b)  experimental  variogram  of  raw 

values,  c)  &  d)  experimental  variogram  and  fitted  models  of  quadratic 
residuals  -  short  wavelength  PEF  and  long  wavelength  PPF  model. 
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a) 


Wavelength  nm 


Wavelength  nm 


b)  c) 


Figure  107.  A.P.  Hill:  pixel  120:  a)  spectrum,  b)  &  c)  experimental  variogram  and 
fitted  models  of  quadratic  residuals  -  short  PES  and  long  wavelength  PPF 
model. 


91 


d) 


e) 


experimental  variogram  and  fitted  models  of  quadratic  residuals  -  short 
wavelength  PEF  and  long  wavelength  PPF,  respectively. 
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a) 


Wavelength  nm 


a) 


Figure  110.  A.P.  Hill:  pixel  141:  a)  spectrum,  b)  &  c)  variogram  and  fitted  models  of 
raw  data,  short  PES  and  long  PPF,  d)  variogram  of  quadratic  residuals  (PPF). 
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a) 


Figure  111.  A.P.  Hill:  pixel  160  a)  spectrum,  b)  experimental  variogram  and  fitted 

model  (PPF)  of  raw  values  c)  experimental  variogram  and  fitted  model 
of  linear  residuals,  long  wavelength  PPF,  and  d)  &  e)  experimental 
variogram  and  fitted  models  of  quadratic  residuals  -  short  wavelength 
PEF  and  long  wavelength  PPF,  respectively. 
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a) 


b)  c) 


Figure  1 12.  A.P.  Hill:  pixel  165  a)  spectrum,  b)  &  c)  experimental  variogram  and 
fitted  models  of  raw  values,  short  and  long  wavelength  d)  &  e) 
experimental  variogram  and  fitted  models  of  linear  residuals  -  short 
wavelength  PEF  and  long  wavelength  PPF,  respectively. 
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a) 


b)  c) 


Figure  113.  A.P.  Hill:  pixel  166:  a)  spectrum,  b)  experimental  variogram  and  fitted 
models  of  raw  data,  short  wavelength  PPF  model  and  c)  experimental 
variogram  and  fitted  models  of  linear  residuals,  short  wavelength  PEF. 


b) 

1 ooooo - 

aoooo- 

60000- 


X 


40000  4 


20000 


*x 


04 - r- 

0  too 


XXxX  X 


c) 


90000 

aoooo 
70000- 
60000 - 
50000- 
40000- 
30000 

20000 4 
10000-1 
0 


\  00  200  300  400  500  600  700 


Figure  114.  A.P.  Hill:  pixel  167:  a)  spectrum,  b)  experimental  variogram  of  raw  data 
c)  experimental  variogram  of  linear  residuals. 
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a) 


Figure  1 15.  A.P.  Hill:  pixel  168:  a)  spectrum,  b)  &  c)  experimental  variogram  and 
fitted  models  (PPF)  of  raw  values,  short  and  long  wavelength  d)  &  e) 
experimental  variogram  and  fitted  models  of  linear  residuals  -  short 
wavelength  PES  and  long  wavelength  PPF,  respectively. 
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Table  8:  Summary  of  trend  analysis  and  model  parameters  of  variograms  for  selected  pixels  for  Fort 
A.P.  Hill 


Pixel 

%  trend  removed  ] 

linear  quadratic'; 

Trend 

Used 

Model  Fitted 

PPF  PES  PEF 

Nugget 

Sill 

Model  Parameters 

Distance  Gradient  Exponent  Wavelength  Amplitude 

Phase 

1 

9.2 

52.3 

Q 

F 

0 

1 .2E+06 

239.6 

248.3 

96874 

2.197 

li 

9.2 

52.3 

Q 

F 

0 

6095 

0.8349 

726.1 

195679 

1.911 

40 

9.9 

48.6 

Q 

F 

0 

69229 

0.4522 

244.1 

114515 

2.229 

II 

9.9 

48.6 

Q 

F 

0 

3552 

0.9490 

783.4 

307081 

1.491 

44 

42.3 

Q 

F 

0 

1.1E+06 

389.8 

250.1 

51297 

2.574 

II 

42.3 

Q 

F 

0 

10679 

0.6948 

558.5 

50849 

IB 

50 

21.4 

40.6 

Q 

F 

0 

1.0E+06 

313.9 

1.2630 

249.6 

46023 

2.620 

II 

21.4 

40.6 

Q 

F 

0 

6392 

0.7865 

693.7 

69732 

2.819 

54 

43.1 

Q 

F 

0 

1 .2E+06 

344.2 

1.1720 

245.0 

52823 

2.267 

II 

43.1 

Q 

F 

0 

7959 

0.7668 

629.0 

66375 

3.075 

65 

19.8 

41.1 

Q 

F 

0 

9.1E+05 

327.4 

1.1590 

249.4 

42694 

2.562 

II 

19.8 

41.1 

Q 

F 

0 

4399 

0.8243 

744.3 

72047 

2.349 

75 

19.1 

40.8 

Q 

F 

0 

1.8E+06 

905.2 

0.7433 

244.7 

57908 

2.392 

76 

20.8 

40.9 

Q 

F 

0 

1.1E+06 

324.3 

250.2 

49225 

2.599 

II 

20.8 

40.9 

Q 

F 

0 

6387 

0.7918 

700.7 

72066 

2.689 

79 

21.3 

40.2 

Q 

F 

0 

1 .2E+06 

313.3 

1.2550 

249.5 

50889 

2.659 

II 

21.3 

40.2 

Q 

F 

0 

6351 

0.8069 

748.2 

86600 

2.562 

87 

18.9 

41.6 

Q 

F 

0 

1 .7E+06 

308.4 

1.2040 

248.9 

75989 

2.538 

II 

18.9 

41.6 

Q 

F 

0 

8137 

0.8317 

750.0 

152570 

2.355 

90 

20.5 

42.0 

Q 

F 

0 

1.1E+06 

326.4 

1.1530 

247.7 

2.453 

II 

20.5 

42.0 

Q 

F 

0 

4657 

0.8408 

750.8 

91291 

2.287 

104 

16.4 

42.9 

Q 

F 

0 

2.0E+06 

301.0 

1.1680 

248.2 

94916 

2.493 

II 

16.4 

42.9 

Q 

F 

0 

10754 

0.8120 

733.6 

181498 

2.347 

107 

19.8 

40.8 

Q 

F 

0 

1.2E+06 

308.6 

1.2530 

248.8 

54333 

2.624 

II 

19.8 

40.8 

Q 

F 

0 

5839 

0.8323 

789.1 

106535 

2.366 

115 

19.3 

43.6 

Q 

F 

0 

2.6E+06 

372.5 

250.2 

119159 

2.619 

II 

19.3 

43.6 

Q 

F 

0 

30773 

0.6661 

521.9 

129584 

2.921 

119 

16.7 

42.1 

Q 

F 

0 

2.7E+06 

295.3 

1.2160 

248.5 

121633 

2.558 

II 

16.7 

42.1 

Q 

F 

0 

14801 

0.8100 

743.8 

242183 

2.400 

120 

41.9 

Q 

F 

0 

2.8E+06 

295.9 

1.1990 

247.4 

126810 

2.493 

II 

41.9 

Q 

F 

0 

15870 

0.8081 

725.4 

258755 

2.421 

133  ... 

17.9 

41 .8 

Q 

F 

0 

2.4E+06 

420.5 

239.7 

89353 

" 

17.9 

41.8 

Q 

F 

0 

16544 

0.7444 

613.2 

122421 

139 

12.6 

46.9 

Q 

F 

0 

2.2E+06 

258.2 

243.8 

146716 

2.075 

" 

12.6 

46.9 

Q 

F 

0 

10405 

0.8425 

714.0 

340547 

1.993 

141 

12.1 

44.8 

Q 

F 

0 

15438 

0.7878 

698.8 

336947 

2.089 

II 

12.1 

44.8 

None 

F 

0 

3.3E+06 

287.7 

1.3160 

259.6 

144847 

3.011 

'• 

12.1 

44.8 

None 

F 

0 

29940 

0.7338 

615.5 

232984 

3.029 

160 

59.1 

58.9 

Q 

F 

0 

1.2E+04 

53.6 

275.7 

1983 

1.894 

II 

59.1 

58.9 

Q 

F 

0 

342.4 

0.6043 

757.5 

5253 

1.210 

n 

59.1 

58.9 

Lin 

F 

0 

1.2E+04 

53.9 

275.9 

1989 

1.890 

II 

59.1 

58.9 

Lin 

F 

0 

325.6 

0.6130 

758.2 

5293 

1.203 

■' 

59.1 

58.9 

None 

F 

0 

1618 

0.3777 

269.8 

2826 

2.279 

Q  quadratic  detrend 
Lin  linear  detrend 


F  is  best  fitting  mode! 

PPF  periodic  with  power  function 
PES  periodic  with  stable  function 
PEF  periodic  exponential  function 


L  fitted  model  with  longer  wavelength 
C  fitted  model  with  shorter  lag  distance 


% 


Table  8:  Summary  of  trend  analysis  and  model  parameters  of  variograms  for  selected  pixels  for  Fort 
A.P.  Hill  (continued) 


Pixel 

%  trend  removed 

linear  quadratic 

Trend 

Used 

Model  Fitted 

PPF  PES  PEF 

Nuqqet 

Sill 

Model  Parameters 

Distance  Gradient  Exponent  Wavelength  Amplitude 

Phase 

165 

34.4 

34.1 

Lin 

F 

0 

3.0E+04 

47.4 

265.3 

4921 

2.134 

II 

34.4 

34.1 

Lin 

F 

0 

3.0E+04 

24.2 

480.2 

7768 

3.090 

ll 

34.4 

34.1 

None 

F 

0 

3.0E+04 

47.4 

265.3 

4921 

2.134 

" 

34.4 

34.1 

None 

F 

0 

544.2 

0.6967 

723.8 

16742 

1.371 

166 

63.9 

63.6 

Lin 

F 

0 

1.0E+04 

62.8 

281.7 

1719 

1.576 

ll 

63.9 

63.6 

Lin 

F 

0 

308.7 

0.5969 

757.1 

3905 

1.162 

II 

63.9 

63.6 

None 

F 

0 

752.5 

0.4882 

275.7 

2444 

1.968 

168 

70.7 

70.6 

Lin 

F 

0 

6.7E+03 

69.3 

1.74 

289.5 

1274 

1.385 

" 

70.7 

70.6 

Lin 

F 

0 

235.7 

0.57 

722.6 

2586 

1.272 

ll 

70.7 

70.6 

None 

F 

0 

282.1 

0.60 

274.0 

1911 

2.153 

II 

70.7 

70.6 

None 

F 

0 

2.581 

1.39 

800.0 

4389 

1.167 

Q  quadratic  detrend 
Lin  linear  detrend 


F  is  best  fitting  model 
PPF  periodic  with  power  function 

PES  periodic  with  stable  function  L  fitted  model  with  longer  wavelength 

PEF  periodic  exponential  function  C  fitted  model  with  shorter  lag  distance 
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Table  9.  Eigenvalues  and  eigenvectors  for  Fort  A.  P.  Hill  using  all  spectra 

Component  PCI  PC2  PC3 

Eigenvalue  102.1  17.4  2.1 

Percentage  variance  81.0  13.8  1.6 
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Table  9.  Continued 
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Table  10  Classes  from  the  non-hierarchical  classification  of  Fort  A.  P.  Hill  spectra  using  63  wavebands  and  the  first  three  principal  components 
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Table  10.  Continued 
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APPENDIX  I 


This  is  a  GenStat  program  to  fit  periodic  models  to  the  variogram.  This  was  written  to 
fit  models  to  the  spectra  in  the  preceding  report.  The  following  functions  can  be  fitted: 

a)  simple  periodic,  b)  periodic  with  phase  shift,  c)  periodic  with  power  function, 
d)  periodic  with  exponential  function,  and  e)  periodic  with  stable  exponential 
function. 

A  file  with  variogram  values  is  appended  so  that  the  program  can  be  tested. 


job  'variogram  fitting' 
scalar  al,  top,  bot,  pi 

scalar  rex,  nugg,  wave,  cfun,  sfun,  grad,  alpha 
calc  pi=3 . 14159 


calc  bot=0 
calc  top=100 
calc  top=top* (top.gt .0) 
text  [nval=2]  title 
text  [nval=l]  label 

frame  window=l , 2 , 3 , 4 ;  ylower=0 . 20 , 0 . 12 , 0 . 32 , 0 . 52 ;  \ 
yupper=0 .8,0.47,0.87,0.87;  \ 

xlower=0 .1,0.48,0.10,0.48;  xupper=0 .9,0.93,0.8,0.93 

variate  [nval=201]  dpar2 ,  dparl,  ypow,  yper,  yabs 
variate  [nval=201]  xlag;  values= !  ( 1 .  .  .  2 01 ) 
variate  [nval=201]  yvar 
variate  [nval=2]  xdv,ydv 
calc  rex=400 

calc  xlag= (xlag-1 . 0) *rex/200 
calc  xdv$ [ 1 ] =0 
calc  xdv$[2]=rex 


n 

open  ' testper . gra ' ;  channel=4 
device  4 


f iletype=graphics 


open  ' 2fthd_37n.gam' ;  channel=2 

read  [channel=2]  title 

print  title 

read  [channel=2]  label 

print  label 

read  [channel=2;  setnvalues=y ;  skip=*]  lag,  gam,  wt 
graph[nr=19;  nc=45]  gam; lag;meth=p 

print  'Model:  unbounded  linear’ 

model  [weights=wt]  gam;  residuals=R;  f ittedvalues=F 
fit  lag 


print  '*****  Periodic  model  with  phase 
print  '  weights  =  npairs  ' 
expression  periodic [ 1 , 2 ] ;  \ 

value=!e(cl  =  ( cos ( 2 *pi*lag/al ) ) } ,  \ 
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DRAFT 


!e(c 2  =  ( sin ( 2*pi*lag/al ) ) ) 

model  [weights^wt;  distribution=normal]  gam;  f ittedvalues=F 
rcycle  al;  initial  =  300;  upper=500;  lower=50^ 

fitnonlinear  [print=monitoring ; calculation=periodic [ ] ]  cl,c2 

rkeep  estimates=kest 
print  kest 
calc  wave=kest$ [ 1] 
calc  cfun=kest$ [2 ] 
calc  sfun=kest$ [3] 
calc  nugg=kest$ [4] 

calc  hypo=sqrt (sfun*sfun+cfun*cfun) 

calc  phi=arccos (cfun/hypo) 

calc  amp=cfun/cos (phi ) 

print  amp,  phi 

calc  yabs=0 

calc  dpar2  =  2  *pi*xlag/wave 

calc  yper=cfun* (cos (dpar2) ) +sfun* (sin(dpar2) ) 

calc  yvar=yper+nugg 

pen  19;  size=1.0 

axes  window=l;  pen=19;  \ 

xtitle='Lag  distance’;  \ 
ytitle= ' Variance '  ;  \ 

ylower=bot ;  yupper=top;  xlower=0;  xupper=rex;  \ 
xmarks= 1(0,100... 300)  ;  ymarks= ! (0,20. . .100) 

ylabels=! t ( 'O',  '1.0',  ’2. O’. '3. O')  " 

"  ylabels= ! t ( ' 0 '  ,  '0.2'  ,'0.4',  '0.6'  ,'0.8',  '1.0')" 

pen  1,2;  linestyle=l , 0 ; colour=l ;method=monotonic , point ;  \ 

symbols=0,4;  thickness=2 . 0 ;  size=0.7 
dgraph  [window^l ; keywindow=0 ; title= ' Variogram ' ]  \ 

y=yvar,gam;  x=xlag,lag;  pen=l,2 


print  ’*****  Periodic  model  *****' 

print  '  weights  =  npairs 
expression  periodic[l];  \ 
value=!e(cl  =  (cos (2*pi*lag/al) ) ) 

model  [weights=wt;  distribution=normal ]  gam;  f ittedvalues=F 

rcycle  al;  initial  =  300;  upper=500;  lower=50 

fitnonlinear  [pr int=monitoring ; calculation=periodic [ ] ]  cl 

rkeep  estimates=kest 

print  kest 

calc  wave=kest$ [ 1 ] 

calc  cfun=kest$ [ 2 ] 

calc  sfun=0 

calc  nugg=kest$ [ 3 ] 

calc  hypo=sqrt (sfun*sfun+cfun*cfun) 

calc  phi=arccos ( cfun/hypo) 

calc  amp=cfun/cos (phi ) 

print  amp,  phi 

calc  yabs=0 

calc  dpar2=2 *pi*xlag/wave 

calc  yper=cfun* (cos (dpar2 ) ) +sfun* (sin (dpar2 ) ) 

calc  yvar=yper+nugg 

pen  19;  size=1.0 

axes  window=l;  pen=19;  \ 

xtitle='Lag  distance';  \ 
ytitle= ' Variance ' ;  \ 

ylower=bot;  yupper=top;  xlower=0;  xupper=rex;  \ 
xmarks=!  (0, 100  .  .  .300)  ;  ymarks= 1  (0 , 20 .  .  . 100)  ;  \ 

"  ylabels= !t('0.2','0.4','0.6','0.8','1.0')" 

pen  1,2;  linestyle=l , 0 ; colour=l ;method=monotonic , point ;  \ 

symbols=0,4;  thickness=2 . 0 ;  size=0.7 


DRAFT 


dgraph  [window=l ; keywindow=0 ; title= ' Variogram' ]  \ 

y =:yvar  ,  gam  ;  x — xlag,lag,  pan  1,2 


print  '***  Periodic  model  with  phase  shift  and  power  function 
print  '  weights  =  npairs  ' 
expression  periodicpower [1 , 2 , 3 ] ;  \ 

value= ! e (cl  =  (cos (2*pi*lag/al) ) ) ,  \ 

!e(c2  =  (sin(2*pi*lag/al)  )  )  ,  \ 

! e (c3= ( ( lag*  *a2 ) ) ) 

model  [weights=wt ;  distribution=normal ]  gam;  f ittedvalues=F 
rcycle  al,a2;  initial  =  400,  1.0;  upper=800 , 1 . 9 ;  lower=5,0.1 
fitnonlinear  [constant=omit; \ 

calculation=periodicpower [ ] ]  cl , c2 , c3 

rkeep  estimates=kest 
print  kest 
calc  wave=kest$ [1] 
calc  cfun=kest$ [ 3 ] 
calc  sfun=kest$ ( 4 ] 
calc  nugg=0 
calc  grad=kest$ [ 5 ] 
calc  alpha=kest$ [ 2 ] 

calc  hypo=sqrt (sfun*sfun+cfun*cfun) 
calc  phi=arccos (cfun/hypo) 
calc  amp=cfun/cos (phi) 

print  amp,  phi,  nugg,  wave,  grad,  alpha 
calc  yabs=0 

calc  dpar2=2*pi*xlag/wave 

calc  yper=cfun* (cos (dpar2 ) ) +sfun* (sin (dpar2 ) ) 
calc  ypow=grad*xlag* *alpha 
calc  yvar=nugg+ypow+yper 

pen  19;  size=1.0 

axes  window=l;  pen=19;  \ 

xtitle='Lag  distance';  \ 
ytitle= 'Variance ' ;  \ 

ylower=bot;  yupper=top;  xlower=0;  xupper=rex;  \ 
xmarks= !  ( 0 , 100  .  .  .  700 )  ;  ymarks= ! ( 0 , 2000 . . . 8 000 ) 

M 

;  \ 

ylabels= !t( ' O' ,  '1.0' ,  '2.0' ,  '3.0' ) 

ii 

pen  1,2;  linestyle=l , 0 ; colour=l ; method=monotonic , point ;  \ 
symbols=0,4;  thickness=2 . 0 ;  size=0.7 
dgraph  [window=l ; keywindow=0 ; titles ' Variogram' ]  \ 

y=yvar,gam;  x=xlag,lag;  pen=l,2 


print  '*****  Periodic  model  with  phase  and  exponential 

function***** 1 

print  ’  weights  =  npairs  1 

expression  periodicexp [ 1 , 2 , 3 ] ;  \ 

value=!e(cl  =  (cos (2*pi*lag/al) ) ) ,  \ 

!e(c2  =  (sin(2*pi*lag/al) ) ) ,  \ 

! e (c3= ( 1 . 0-exp ( - ( lag/a2 ) ) ) ) 

model  [weights=wt;  distribution=normal]  gam;  f ittedvalues=F 
rcycle  al,a2;  initial  =  400,100;  upper=900 , 500 ;  lower=5,l 


DRAFT 


fitnonlinear  [ constant=omit ; calculation=periodicexp [ ] ]  \ 

cl , c2 , c3 

rkeep  estimates=kest 
print  kest 
calc  wave=kest$ [ 1 ] 
calc  cf un=kest$  [ 3 ] 
calc  sfun=kest$ [4] 
calc  nugg=0 
calc  dist=kest$  [2] 
calc  sill=kest$ [5] 


calc  hypo=sqrt (sfun*sfun+cfun*cfun) 
calc  phi=arccos (cfun/hypo) 
calc  amp=cfun/cos (phi ) 

print  sill,  dist,  amp,  cfun,  sfun,  wave 
calc  yabs=0 

calc  dpar2  =  2  *pi*xlag/wave 

calc  yper=cfun* (cos (dpar2) ) +sfun* (sin(dpar2) ) 

calc  dparl=xlag/dist 

calc  yexp=sill* (1 . 0-exp (-dparl) ) 

calc  yvar=nugg+yexp+yper 

pen  19 ;  size=l . 0 

axes  window=l;  pen=19;  \ 

xtitle='Lag  distance';  \ 
ytitle= 1  Variance '  ;  \ 

ylower=bot ;  yupper=top;  xlower=0;  xupper=rex;  \ 
xmarks= ! ( 0 , 100 . . . 700)  ;  ymarks= !  (0 , 1000  ...  7000) 

"  ylabels=!t('0' , '1.0', '2.0' , '3.0')  " 

pen  1,2;  linestyle=l , 0 ; colour=l ;method=monotonic , point ;  \ 

symbols=0,4;  thickness=2 . 0 ;  size=0.7 
dgraph  [window=l ; keywindow=0 ; titles  1 Variogram ' ]  \ 

y=yvar,gam;  x=xlag,lag;  pen=l,2 


print  ’***  Periodic  model  with  phase  shift  and  stable  function 
print  1  weights  =  npairs 
expression  periodicexpst [ 1 , 2 , 3 ] ;  \ 

value=!e(cl  =  ( cos ( 2 *pi * lag/al ) ) ) ,  \ 

!e(c2  =  (sin (2*pi*lag/al) ) )  ,  \ 

! e (c3= ( 1 . 0-exp ( - ( lag/a2 )  )  )  ) 

model  [weights=wt;  distr ibution=normal ]  gam;  f ittedvalues=F 
rcycle  al,a2;  initial  =  400,100;  upper=800 , 1000 ;  lower=5,l 
fitnonlinear  [ constant=omit ;  \ 

calculation=periodicexpst [ ] ]  cl , c2 , c3 

print  '***  Periodic  model  with  phase  shift  and  stable  function 
print  '  weights  =  npairs  ’ 
expression  periodicexpst [ 1 , 2 , 3 ] ;  \ 

value=!e(cl  =  (cos (2*pi*lag/al) ) ) ,  \ 

!e(c2  =  ( sin (2*pi*lag/al ) ) ) ,  \ 

! e (c3= ( 1 . 0-exp ( - ( (lag**a3) / (a2**a3) ) ) ) ) 

model  [weights=wt;  distribution=normal ]  gam;  f ittedvalues^F 
rcycle  al,a2,a3;  initial  =  400,100,1;  upper=800 , 1000 , 2 ; 
lower =5 , 0 . 1 , 0 . 1 ; \ 
step=0 .1,0.1,0.01 
fitnonlinear  [const=omit ;  \ 
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calculation=periodicexpst [ ] ]  cl , c2 , c3 
rkeep  estimates^kest 
print  kest 
calc  wave=kest$ [1] 
calc  cfun=kest$ [4] 
calc  sfun=kest$ [5] 
calc  nugg=0 
calc  alpha=kest$ [3 ] 
calc  dist=kest$ [2] 
calc  sill=kest$ [6] 


calc  hypo=sqrt (sfun*sfun+cfun*cfun) 
calc  phi=arccos (cfun/hypo) 
calc  amp=cfun/cos (phi ) 

print  sill,  dist,  amp,  cfun,  sfun,  'wave,  alpha 
calc  yabs=0 

calc  dpar2=2 *pi*xlag/wave 

calc  yper=cfun* (cos (dpar2) ) +sfun* (sin(dpar2) ) 
calc  dparl=xlag**alpha/dist**alpha 
calc  yexps=sill* (1 . 0-exp (-dparl) ) 

calc  yvar=nugg+yexps+yper 

pen  19;  size=1.0 

axes  window=l;  pen=19;  \ 

xtitle='Lag  distance';  \ 
ytitle= 'Variance ' ;  \ 

ylower=bot;  yupper=top;  xlower=0;  xupper=rex;  \ 
xmarks= !  ( 0 , 100  .  .  .  700)  ;  ymarks= !  ( 0 , 1000  .  .  . 8 000 ) 

"  ylabels= !t('0',  '1.0',  '2.0',  '3.0')  " 

pen  1,2;  linestyle=l , 0 ; colour=l ; method=monotonic , point ;  \ 
symbols=0 , 4 ;  thickness=2 . 0 ;  size=0.7 
dgraph  [window=l ; keywindow=0 ; title= 1 Variogram 1 ]  \ 

y=yvar,gam;  x-xlag,lag;  pen=l,2 


stop 


Variogram  values 


'  Lag[0]  Vgram[0]  Count [0] 

'  Fort  Hood  2 ' 

'  Position  37 ' 

12.0  3.52  234.0 
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APPENDIX  II 

Paper  submitted  to  International  Journal  of  Remote  Sensing  -  please  do  not  quote 
from  this  or  use  the  information  therein. 


Using  the  Variogram  to  Explore  Imagery  of  Two  Different  Spatial  Resolutions 

M.  A.  OLIVER1,  J.  A.  SHINE2  and  K.  R.  SLOCUM2 

1  Department  of  Soil  Science,  University  of  Reading,  Whiteknights,  Reading  RG6 
6DW,  United  Kingdom. 

2United  States  Army  Topographic  Engineering  Center,  7701  Telegraph  Road, 
Alexandria,  Virginia  22315-3864,  USA 

Abstract.  The  resolution  of  remotely  sensed  data  is  becoming  increasingly  fine,  and 
there  are  now  many  sources  of  data  with  a  pixel  size  of  lmxlm.  This  produces 
huge  amounts  of  data  that  have  to  be  stored,  processed  and  transmitted.  For 
environmental  applications  this  resolution  possibly  provides  far  more  data  than  are 
needed;  data  overload.  This  poses  the  question:  how  much  is  too  much?  We  have 
explored  two  resolutions  of  data,  1-m  pixel  CAMIS  data  and  20-m  pixel  SPOT  data, 
using  the  variogram  of  geostatistics.  For  both  we  used  the  normalized  difference 
vegetation  index  (NDVI).  It  is  clear  that  there  is  more  information  in  the  1-m  data  and 
three  scales  of  spatial  variation  were  identified:  a  short  scale  of  about  25  m,  an 
intermediate  one  of  about  150  m  and  a  much  longer  one  of  about  600  m.  We  sub¬ 
sampled  the  1-m  data  and  identified  an  additional  scale  of  variation  of  about  300  m. 
The  latter  and  the  short  spatial  scale  were  identified  consistently  until  the  sub-sample 
was  1  pixel  in  15  for  each  row  and  column  (or  1  pixel  in  225).  At  this  stage  spatial 
scales  of  about  100  m  and  600  m  were  described,  which  suggested  that  only  now  was 
there  a  real  difference  in  the  amount  of  information  available.  We  compared  the 
results  of  this  analysis  with  the  variogram  of  the  SPOT  data.  Similar  spatial  scales  to 
those  from  the  sub-sampled  1-m  data  were  identified.  From  this  analysis  it  seems  that 
a  pixel  size  of  20-m  is  adequate  for  many  environmental  applications,  and  that  if  more 
detail  is  required  the  1-m  data  could  be  sub-sampled  to  1  pixel  in  10  for  every  row  and 
column  without  any  serious  loss  of  information.  This  reduces  significantly  he  amount 
of  data  the  needs  to  be  stored,  transmitted  and  analysed  and  has  important 
implications  for  data  compression. 


1.  Introduction 

The  applications  of  remotely  sensed  imagery  are  constantly  increasing  as  are  the 
capabilities  for  increasingly  fine  spatial  resolution.  Aerial  images  contain  large 
amounts  of  information,  both  spatial  (pixels)  and  spectral  (wavebands)  data.  As  the 
resolution  increases  this  results  in  vastly  more  information:  imagery  with  a  spatial 
resolution  of  1  m  x  1  m  per  pixel  contains  400  times  as  much  data  as  that  with  a  pixel 
size  of  20  m  x  20  m.  A  wide  variety  of  satellite  sensors  will  provide  information  at 
1-m  or  near  1-m  resolution  in  the  next  few  years  (Aplin  et  al.,  1997).  Sensor  systems 
such  IKONOS  (Space  Imaging,  Corp.,  panchromatic  channel)  and  CAMIS 
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(Computerized  Airborne  Multispectral  Imaging  System)  already  achieve  this 
resolution  (Carder,  1993;  Baker,  1999).  Hyperspectral  imagery,  such  as  HYDICE.  has 
many  more  wavebands  than  the  conventional  imagery  of  SPOT  (4  bands)  or  Landsat 
TM  (7  Bands).  For  example  SEBASS  thermal  hyperspectral  imagery  has  128 
wavebands.  With  such  data,  if  the  spatial  resolution  is  the  same  as  that  for  a  SPOT 
image,  then  there  are  almost  42  times  more  data.  In  this  paper  we  considei  the  spatial 
resolution  of  the  data  only. 


If  both  the  resolution  and  spectral  information  are  increased  then  it  is  clear  how  the 
amount  of  information  will  expand.  There  are  consequences  of  this  for  data  storage, 
transfer  and  analysis.  It  is  evident  that  1-m  imagery  contains  more  information  than 
that  at  coarser  resolutions,  such  as  the  20-m  SPOT  data,  but  this  raises  two  questions. 
How  much  of  the  additional  information  is  of  real  use?  What  are  the  relative  costs  of 
obtaining,  analysing  and  storing  the  extra  data?  Acquiring  and  storing  1-m  data  may 
entail  a  significant  cost,  even  in  an  era  of  cheaper  hardware  and  constantly  increasing 
memory  availability.  The  speed  and  accuracy  with  which  data  can  be  transferred  and 
processed  will  be  improved  by  reducing  the  data  optimally  for  future  use.  The  extent 
of  data  compression  will  depend  on  the  degree  of  redundancy  present  in  the  data  and 
what  the  future  uses  are  likely  to  be.  By  examining  1-m  imagery  at  a  range  of 
resolutions  it  might  determine  the  stages  at  which  crucial  spatial  information  is  lost 
and,  therefore,  whether  processing  the  full  set  of  data  is  worthwhile. 


Geostatistics  has  been  applied  successfully  to  remotely  sensed  data  in  a  range  ways, 
for  example  to  improve  image  classification  (Abarca-Hernandez  and  Chica-Olmo, 
1999),  to  monitor  crop  growth  (Erickson  et  al,  1999),  and  to  design  optimal  sampling 
schemes  for  ground  surveys  (Oliver  et  al.,  2000),  In  this  paper  we  use  the  variogram 
of  Regionalized  Variable  Theory  to  compare  the  data  from  a  SPOT  satellite  image  and 
1-m  CAMIS  imagery  over  eastern  Virginia  in  the  United  States.  The  1-m  data  were 
sub-sampled  to  explore  a  range  of  distances  between  pixels,  but  less  than  that  of  the 
SPOT  20-m  pixel.  The  analysis  used  the  normalized  difference  vegetation  index 
(NDVI)  for  both  sets  of  data.  The  NDVI  is  related  to  the  proportion  of 
photosynthetically  absorbed  radiation  and  it  is  therefore  useful  for  identifying  changes 
in  ground  cover  and  possibly  separating  landscape  features.  Chavez  (1992)  showed 
that  the  NIR  waveband  contains  greater  within-scene  variation  than  the  red  green  or 
blue  wavebands.  Since  NDVI  incorporates  NIR  maximum  variability  is  ensured. 


2.  The  study  area  and  image  data 

The  area  studied  is  part  of  Fort  A.  P.  Hill  in  northeastern  Virginia.  It  is  located  on  the 
gently  undulating  terrain  of  the  Atlantic  Coastal  Plain  of  the  United  States  of  America. 
The  site  was  selected  because  of  its  complex  topography  which  gives  rise,  in  part,  to 
the  diversity  in  land  cover  in  a  small  geographic  area.  The  study  area  has  been 
intensely  dissected  by  many  small  ephemeral  drainage  channels  that  help  to  recharge 
the  many  ponds  and  wetlands,  but  the  average  relief  is  only  about  15  m.  Some  of  these 
ponds  are  evident  as  blue  areas  in  the  pixel  map  of  NDVI  from  the  SPOT  image, 
Figure  1.  Drainage  is  to  two  main  river  basins:  the  Rappahannock  and  the  James. 
Hence  there  is  a  significant  watershed  crossing  the  Fort  in  an  W-E  direction.  The  area 
has  mixed  land  use,  the  most  extensive  of  which  are  the  hardwood,  pine  and  mixed 
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forests.  In  addition  there  is  maintained  grassland,  shrub,  and  man-made  features  such 
as  buildings  and  roads. 

2.1.  The  images 

The  image  obtained  from  the  SPOT  satellite  has  three  wavebands,  green  (G),  red  (R) 
and  near  infrared  (NIR).  The  normalized  difference  vegetation  index  was  calculated 
from  the  red  and  NIR  wavebands  using  the  usual  equation.  The  part  of  the  SPOT 
image  analysed  had  178  rows  and  151  columns,  giving  a  total  of  26  878  pixels.  It 
covered  an  area  of  3.6  km  x  3.0  km.  The  20-m  pixel  size  of  the  SPOT  image  was  the 
largest  spatial  resolution  examined  in  this  analysis.  The  1  m  x  1  m  pixel  image 
obtained  from  CAMIS  sensors  has  four  wavebands,  blue  (B),  green  (G),  red  (R)  and 
NIR.  As  for  SPOT  the  NDVI  was  calculated.  There  were  some  gaps  in  the  cover  of 
this  image  because  it  was  formed  from  a  photo-mosaic,  comprising  8  499  498  pixels, 
Figure  2.  Both  images  covered  more  or  less  the  same  area.  To  explore  fully  the  range 
of  spatial  structures  that  might  be  present  at  different  scales,  the  1-m  data  was  sub¬ 
sampled  to  provide  information  at  an  increasingly  large  separation.  At  the  same  time 
this  reduced  the  amount  of  data  for  analysis  substantially.  Another  way  of  achieving 
this  would  have  been  to  average  the  digital  information  over  the  relevant  block  of 
pixels.  This  would  have  changed  the  support  of  the  original  data  and  would  have 
reduced  the  variance  present.  Since  our  aim  was  to  examine  the  extent  of  redundancy 
in  the  data  we  did  not  modify  the  original  data  by  increasing  the  pixel  size.  The  effect 
of  this  is  something  that  should  be  explored  in  the  future. 

Figure  2  shows  the  effect  of  combining  several  separate  1-m  images  to  cover  the  study 
area;  the  colour  intensity  is  not  uniform.  Since  this  could  affect  the  variogram  analysis 
we  limited  the  sub-sampling  to  a  maximum  separation  of  1  pixel  in  every  15  along 
each  row  and  column  (or  1  pixel  in  a  block  of  225).  This  ensured  that  the  comparisons 
between  pixels  did  not  extend  over  distances  much  larger  than  1000  m  (about  a  third 
of  the  extent  of  the  area). 

The  first  sub-sample  retained  one  pixel  in  every  two  along  each  row  and  column.  This 
gave  a  2-m  separation  between  the  centroids  of  neighbouring  pixels.  This  achieved  a 
75  %  reduction  in  the  data.  This  process  was  repeated  to  give  a  3-m  separation 
between  pixels  (89  %.  reduction  in  data),  a  6-m  separation  (97.22  %  reduction  in 
data),  pixels  10-m  apart  (99  %  reduction  in  data),  and  pixels  15-m  apart  (99.6  % 
reduction  in  data).  After  the  last  reduction  there  were  37  553  pixels  in  the  data,  which 
was  still  more  than  the  number  of  SPOT  pixels. 


3  Methodology 

Spatial  properties  vary  from  place  to  place  and  this  can  occur  at  more  than  one  scale 
simultaneously.  As  a  result  the  spatial  structures  can  vary  at  scales  that  differ  by 
several  orders  of  magnitude,  from  a  few  metres  to  hundreds  of  kilometres.  This  is 
known  as  nested  variation.  Such  variation  in  spatial  properties,  including  that  in 
spectral  imagery,  can  be  described  using  the  theory  and  methods  of  geostatistics 
(Oliver  et  al,  2001).  The  variogram  is  the  central  tool  of  this  methodology;  it 
provides  an  unbiased  description  of  the  scale  and  pattern  of  spatial  variation.  It  does 
so  by  measuring  the  degree  of  spatial  correlation  between  sampling  points,  based  on 
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the  notion  that  the  property  is  more  likely  to  be  similar  at  nearby  locations  than  at 
distant  ones  (Webster  and  Oliver,  2001).  The  variogram  is  defined  by: 


r(h)=^E[{z(x)-Z(x  +  h)}2]  0) 

where  Z(x)  and  Z(x  +  h)  are  the  values  of  Z,  the  property  of  interest,  at  any  two  places 
x  and  x  +  h  separated  by  h,  which  is  a  vector  in  both  distance  and  direction,  known  as 
the  lag.  The  symbol  E  denotes  the  expectation.  The  semivariance,  y  at  a  given 
separation  is  half  the  expected  squared  difference  between  values  at  that  separation.  It 
summarises  the  way  in  which  the  variance  of  a  property  changes  as  the  distance  and 
direction  separating  any  two  points  varies. 


We  used  variography  to  explore  the  spatial  information  in  the  SPOT  and  1-m  image 
data.  The  standard  computing  formula  for  estimating  the  semivariance  at  a  given  lag, 
h  is: 


?(h)=-^rrf M*.  +»»2 


2M  (h)  £ 


(2) 


where  y  (h)  is  the  estimate  of  7(h),  z(x\)  and  z(x j  +  h)  are  the  observed  values  of  Z  (the 
digital  numbers  in  this  case)  at  x,  and  Xj  +  h  respectively,  and  M(h)  is  the  number  of 
paired  comparisons  at  that  lag.  By  changing  h  an  ordered  set  of  values  is  obtained,  and 
this  is  the  sample  or  experimental  variogram. 

Mathematical  functions  or  models  were  fitted  to  the  experimental  variogram  by 
weighted  least  squares  approximation  using  GenStat  (Payne,  2000).  The  parameters  of 
the  variogram  model  describe  the  structure  of  the  spatial  variation  and  the  scale(s)  of 
variation  present. 


4.  Analysis  and  results 

Experimental  variograms  were  computed  from  the  digital  numbers  of  the  three 
wavebands  of  the  SPOT  image  and  for  the  normalized  difference  vegetation  index 
(NDVI).  Since  the  results  for  all  are  similar  we  describe  those  for  NDVI  only. 
Variograms  were  computed  along  the  rows  and  columns  of  the  SPOT  image,  and 
these  were  then  averaged  to  give  the  overall  variogram.  They  were  computed  to  a 
maximum  lag  of  100  pixels  (2000  m)  for  SPOT  imagery  and  1000  pixels  (1000  m)  for 
CAMIS  imagery. 


Figure  3a  shows  the  average  experimental  variogram  for  NDVI  from  the  SPOT  image 
computed  to  a  lag  of  100.  It  is  evident  that  after  a  lag  of  about  40  pixels  the 
experimental  semivariances  become  erratic  and  after  a  lag  of  80  pixels  they  start  to 
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decrease.  Since  both  of  these  features  can  have  an  adverse  effect  on  modelling  the 
variogram  we  fitted  models  to  maximum  lags  of  40  and  80  pixels.  A  range  of 
authorized  mathematical  functions  was  tried;  this  range  is  somewhat  limited  because 
they  must  be  conditional  negative  semi-definite  (see  Webster  and  Oliver,  2001).  The 
best  fitting  models  to  the  experimental  variograms,  in  the  least  squares  sense,  were 
nested  exponential  and  spherical  functions,  which  describe  more  than  one  spatial 
scale.  Double  exponential  and  spherical  functions  fitted  all  experimental  variograms, 
except  for  that  for  the  1-m  lag,  which  was  fitted  best  by  a  triple  spherical  function. 
The  equation  for  the  double  spherical  function  is: 


y(h)  =  c0  +  Cj 


3 h 

2a  i  2 


yUl  j 


+  C , 


3  h  1 


2  a2  2 


r_h  A 

Kai  J 


r{h)=c0+cl+cr 
y{h)  =  c0+cl+c2 


3  h  1 


2  a2  2 


(  u  \ 
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\  z  J 


for  h<al 


for  h<a2 
for  h>  a2 

(3) 


where  c0is  the  nugget  variance,  c\  is  the  autocorrelated  variance  or  sill  and  ax  is  the 
correlation  range  or  range  of  spatial  dependence  of  the  first  spatial  structure,  and  C2  is 
the  autocorrelated  variance  and  a-i  is  the  range  of  the  second  spatial  structure. 


The  triple  spherical  function  is  given  by: 

for  h<a1 

for  h<a2 

for  h<a3 

for  h>  a3 

(4) 

where  C3  is  autocorrelated  variance  and  a 3  is  the  range  of  the  third  spatial  structure. 

The  nested  exponential  function  is  given  by: 

y{h)=c0  +  c1{l-exp(-h/7j)}  +  c2{l-exp(-h/r2)}  (5) 

where  r\  and  r2  are  the  distance  parameters  of  the  short-range  and  long-range  spatial 
components,  respectively.  The  exponential  model  reaches  its  sill  asymptotically  and 
therefore  does  not  have  a  definite  range  of  spatial  dependence;  a  working  range  can  be 
obtained  as  a  -  3r. 
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In  the  first  instance  we  fitted  double  spherical  and  exponential  functions  to  the 
average  experimental  variogram  computed  to  a  maximum  lag  of  80  pixels.  The  best 
fitting  model  was  the  nested  exponential  one,  which  has  an  asymptotic  sill  as  shown 
in  Figures  3b:  the  symbols  are  the  experimental  values  and  the  solid  line  is  the  fitted 
model.  The  double  spherical  model  is  not  such  a  good  fit  as  shown  in  Figure  3f  and  by 
the  residual  mean  square  (RMS)  in  Table  2.  The  parameters  of  the  fitted  models  are 
given  in  Table  2. 

Fitting  models  to  variograms  is  one  of  the  most  controversial  aspects  of  geostatistics; 
there  are  no  hard  and  fast  rules  and  no  consensus  on  whether  to  fit  by  eye  or 
statistically.  We  prefer  the  latter  because  it  provides  a  means  of  objective  comparison 
between  different  models.  From  experience  we  know  that  model  fitting  is  affected  by 
any  erratic  behaviour  in  the  experimental  values.  Figures  3a  and  b  show  that  after  a 
lag  of  about  40  pixels  the  variograms  fluctuate  about  their  sills. 

On  closer  inspection  of  the  variogram  computed  to  a  lag  80  pixels  it  seemed  more 
sensible  to  model  the  experimental  variogram  to  a  maximum  lag  of  40  pixels  to 
achieve  a  better  fit.  Table  1  gives  the  best  fitting  model  to  the  rows,  columns  and  the 
average  variograms  to  lag  40  pixels,  Figure  3  c,  d,  and  e,  respectively.  A  double 
spherical  provided  the  best  fit  in  all  cases  because  these  variograms  flatten  to  a 
distinct  sill.  Table  2  gives  the  model  parameters;  particular  attention  should  be  paid  to 
the  short-  and  long-range  spatial  structures.  The  variogram  of  the  columns  suggests 
that  both  spatial  components  have  shorter  ranges  in  an  E-W  direction  than  those  for 
the  rows  in  the  N-S  direction.  Thus  the  variation  is  anisotropic.  However,  Figure  1  of 
the  NDVI  values  does  not  suggest  that  this  is  pronounced.  Four  spatial  scales  of 
variation  emerge  from  these  variograms:  about  130  m,  340  m,  450  m,  and  640  m. 

We  analysed  all  of  the  wavebands  of  the  1-m  imagery,  but  we  give  the  results  for 
NDVI  only.  We  also  computed  variograms  of  the  rows,  columns  and  average,  but  we 
present  the  results  of  the  average  variograms  only.  The  observed  differences  between 
those  of  the  rows  and  columns  were  similar  to  those  for  the  SPOT  data.  We  computed 
the  variogram  from  the  1-m  data  to  a  maximum  lag  of  1000  m,  but  fitted  the  model  to 
800  m,  Figure  4a.  The  double  functions  fitted  to  the  SPOT  data  did  not  provide  such  a 
god  fit  and  triple  functions  were  also  fitted.  The  triple  spherical  model  provided  the 
best  fit,  Figure  4a  shows  the  experimental  variogram  and  best  fitting  model,  and  Table 
2  gives  the  model  parameters.  The  double  exponential  model  provided  the  best  fit  to 
all  of  the  variograms  of  the  subsets. 

The  triple  spherical  variogram  model  fitted  to  the  pixels  separated  by  1-m  has  a  short- 
range  component  that  matches  that  of  the  variograms  from  the  first  five  sub-samples 
reasonably  (see  Table  3).  The  middle-range  component  matches  that  of  the  short- 
range  component  identified  in  the  SPOT  image,  and  the  longest  range  matches  that  of 
the  longest  range  fitted  to  both  the  SPOT  data  and  to  the  1  in  15  sub-sample.  The 
variograms  of  the  sub-samples  are  remarkably  consistent  except  for  the  largest  pixel 
spacing  of  the  1  in  15  sub-sample  (1  pixel  in  a  block  of  225),  Figure  4g.  The  scale  of 
spatial  dependence  of  the  short-range  component  for  the  2-m,  3-m,  4-m,  6-m  and  10- 
m,  Figure  4b,  c,  d,  e  and  f  separations  is  about  33  m,  and  that  for  the  long-range  one  is 
between  310  m  and  346  m.  There  is  no  equivalent  of  this  short-range  component  in 
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the  SPOT  data  because  it  is  close  to  the  size  of  the  pixels.  For  the  longer-range 
component  this  is  close  to  the  short-range  parameter  of  the  average  variogram 
computed  to  the  longer  maximum  lag  (80  pixels),  Table  2.  For  the  largest  separation 
between  pixels  of  1  in  15  the  long-  and  short-range  components  are  similar  to  those  of 
the  average  variogram  for  the  shorter  maximum  lag  (40  pixels)  for  the  SPOT  data  and 
to  the  largest  range  fitted  to  the  1-m  data. 

6.  Discussion  and  Conclusions 

For  these  particular  data  one  outstanding  conclusion  from  the  results  is  that  for  the 
1 -metre  CAMIS  imagery  there  is  little  to  be  gained  in  terms  of  detecting  spatial 
structures  of  environmental  interest  with  pixels  less  than  10-m  apart.  Therefore,  the 
original  1-m  data  could  be  compressed  by  1: 100  in  2-D  space,  reducing  the  data  by  m 
99  %  without  affecting  the  variogram  results.  This  suggests  that  substantial 
compression  of  huge  1 -metre  imagery  data  sets  is  possible  without  losing  significant 
spatial  information.  On  the  other  hand  if  information  about  man-made  features  or 
individual  trees  is  of  interest  then  this  will  not  be  so,  but  for  imagery  that  is  to  be  used 
for  monitoring  the  landscape  and  other  environmental  tasks,  this  degree  of  data 
compression  would  be  feasible. 

The  results  of  the  1-m  data  analysis  are  remarkably  consistent  in  terms  of  the  spatial 
scales  that  they  identify,  and  the  longer  ones  match  with  scales  of  variation  identified 
from  the  SPOT  data.  Nevertheless,  the  analysis  of  the  1-m  imagery  identifies  a  scale 
of  spatial  variation  that  is  not  evident  from  the  SPOT  data  of  the  same  region  (Fort 
A.P.  Hill,  Virginia,  USA).  This  scale  of  variation  of  about  30-m  appears  because  of 
the  increased  resolution  of  the  CAMIS  imagery.  It  can  be  seen  on  the  imagery,  Figure 
2,  and  appears  as  ‘pock  marks’.  It  seems  likely  that  this  is  a  shadowing  effect  from  the 
trees.  The  scales  of  variation  identified  in  the  sub-sample  of  1  pixel  in  every  15  rows 
and  every  15  columns  are  similar  to  those  identified  from  the  SPOT  data  when  the 
maximum  lag  to  which  the  variogram  was  calculated  is  40  pixels.  This  finding  is 
significant  as  researchers  attempt  to  move  seamlessly  from  one  scale  of  spatial 
resolution  to  another  to  satisfy  the  needs  in  multi-purpose  landscape  characterization. 

For  the  SPOT  data  modelling  the  variogram  was  more  difficult  because  of  fluctuations 
in  the  experimental  semivariances  once  the  sill  had  been  reached.  Several 
permutations  of  modelling  these  variograms  were  tried  and  the  results  differed  more 
or  less  (see  Figure  3b  and  f).  The  approach  to  modelling  can  vary  according  to  what 
the  model  is  to  be  used  for.  A  principal  reason  for  modelling  the  variogram  is  to 
provide  the  spatial  information  needed  for  geostatistical  prediction,  kriging.  In  this 
case  the  fit  of  the  variogram  close  to  the  origin  is  the  most  important  because  the  first 
few  lags  only  are  generally  used  in  local  estimation.  For  exploring  the  spatial  scales  of 
variation  present  the  overall  fit  of  the  model  to  different  maximum  lags  becomes 
important  and  that  is  the  approach  that  we  have  adopted  here. 

A  geostatistical  comparison  between  information  at  different  pixel  separations 
provides  an  opportunity  for  deciding  what  level  of  data  compression  would  be 
acceptable  in  an  objective  way.  This  could  be  investigated  further  by  examining  the 
quality  of  predictions  from  different  methods  of  data  restoration  such  as  kriging  and 
wavelet  analysis  using  a  cross  validation  procedure  (Oliver,  et  al.,  2000). 
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Figure  Captions 

Figure  1.  Pixel  map  of  NDVI  from  a  SPOT  image  for  part  of  Fort  A.  P.  Hill,  Virginia, 
USA. 

Figure  2.  Pixel  map  of  1-m  CAMIS  imagery  for  part  of  Fort  A.  P.  Hill,  Virginia, 

USA. 

Figure  3.  Experimental  variograms  (symbols)  and  fitted  models  (solid  lines)  of  NDVI 
from  the  SPOT  data:  a)  average  of  rows  and  columns  computed  to  100  lags  (2000  m), 
b)  average  of  rows  and  columns  computed  to  80  lags  fitted  by  an  exponential  function 
(1600  m),  c)  rows  computed  to  40  lags  (800  m),  d)  columns  computed  to  40  lags  (800 
m),  e)  average  computed  to  40  lags  (800  m),  f)  average  of  rows  and  columns 
computed  to  80  lags  fitted  by  a  spherical  function  (1600  m). 

Figure  4.  Experimental  variograms  (symbols)  and  fitted  models  (solid  lines)  of  NDVI 
from  the  CAMIS  data:  a)  1-m,  b)  2-m,  c)  3-m,  d)  4-m,  e)  6-m,  f)  10-m  and  g)  15-m 
separation  between  pixels. 
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Table  1.  Some  current  sensor  systems  with  their  spatial  resolution  and  number  of 
bands 


Sensor  System 

Spatial  Resolution 

#  wavebands 

Computerized  Airborne  Multispectral  Imaging  System 

1-m  or  less 

4 

Space  Imaging-IKONOS 

4-mMSI;  1-m  Pan 

4,  1 

SPOT 

20-m  MSI;  10-m  Pan 

3,  1 

Landsat  Thematic  Mapper 

30-m 

7 

HYDICE 

0.8  to  4-m 

210 

Airborne  Visible/Infrared  Imaging  Spectrometer 

20-m 

224 

Advanced  Very  High  Resolution  Radiometer 

1.1-km 

5 

Table  2.  Model  parameters  for  variograms  from  the  SPOT  Imagery 


Variogram 

Model 

RMS 

co 

Cl 

C2 

fli(m) 

a2(m) 

Rows  40 

Double  exponential 

27.61 

0 

0.0341 

0.0560 

322.5 

10602 

Columns  40 

Double  exponential 

No  suitable  fit  achieved  for  this  model 

Average  40 

Double  exponential 

63.21 

0 

0.0360 

0.0080 

344.0 

1400.4 

Average  80 

Double  spherical 

52.39 

0 

0.0214 

0.0225 

139.2 

674.8 

Table  3.  Model  parameters  for  the  average  variograms  of  different  pixel  separations 
for  the  1-m  imagery 


Variance 

Range 

Data  set 

Model  type 

C0 

<4 

C2 

c,  axl m 

ajm 

cz,/m 

1  metre  data 

Triple 

0 

0.01728 

0.02560 

0.01557  23.42 

147.7 

620.8 

1  in  2  pixels 

spherical 

Double 

0 

0.0245 

0.0318 

30.08 

317.1 

1  in  3  pixels 

exponential 

Double 

0 

0.0255 

0.0325 

33.79 

346.2 

1  in  4  pixels 

exponential 

Double 

0 

0.0250 

0.0323 

33.74 

343.3 

1  in  6  pixels 

exponential 

Double 

0 

0.0245 

0.0325 

31.82 

320.2 

1  in  10  pixels 

exponential 

Double 

0 

0.0243 

0.0326 

31.50 

310.2 

1  in  15  pixels 

exponential 

Double 

exponential 

0.0215 

0.0354 

0.0068 

109.1 

626.4 
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APPENDIX  III 

Paper  submitted  to  International  Journal  of  mathematical  Geology  -  please  do  not 
quote  from  this  or  use  the  information  therein. 


Haar  wavelets  and  the  generalization  of  the  variogram 


E.  H.  Bosch1,  M.  A.  Oliver2  and  R.  Webster  3 


The  experimental  variogram  computed  in  the  usual  way  by  the  method  of  moments 
and  the  Haar  wavelet  transform  are  similar  in  that  they  filter  data  and  yield 
informative  summaries  that  may  be  interpreted.  The  variogram  filters  out  the  mean; 
wavelets  can  filter  variation  at  several  spatial  scales  and  thereby  provide  a  richer 
repertoire  for  analysis  and  demand  no  assumptions  other  than  that  of  finite  variance. 
This  paper  compares  the  two  functions,  identifying  that  part  of  the  Haar  wavelet 
transform  that  gives  it  its  advantages,  and  suggests  how  the  variogram  might  be 
elaborated  with  additional  filter  coefficients  to  reveal  features  of  the  data  that  are  not 
evident  in  its  usual  form. 

Two  examples  in  which  soil  data  recorded  at  regular  intervals  on  transects  are 
analysed  illustrate  the  extended  form  of  the  variogram.  The  apparent  periodicity  of 
gilgais  in  Australia  seems  to  be  accentuated  as  filter  coefficients  are  added,  but 
otherwise  provides  no  new  insight.  Adding  filter  coefficients  in  the  analysis  of  the 
topsoil  across  the  Jurassic  scarplands  of  England  changed  the  upper  bound  of  the 
variogram;  it  then  resembled  the  within-class  variogram  computed  by  the  method  of 
moments.  To  elucidate  these  results  we  simulated  several  series  of  data  to  represent 
data  with  no  stepped  transitions,  data  with  stepped  transitions,  data  with  long-range 
linear  trend,  and  data  with  local  trend.  The  results  suggest  that  the  wavelet  variogram 
can  filter  out  the  effects  of  transitions  from  one  class  to  another,  as  across  boundaries, 
and  long-range  trend,  but  not  local  trend. 


KEY  WORDS:  filtering,  Haar  wavelet,  variogram,  wavelet  analysis. 


1  USACE-ERDC-TEC,  7701  Telegraph  Road,  Alexandria,  VA  22315-3864,  USA 

2  Department  of  Soil  Science,  The  University  of  Reading,  P  O  Box  233,  Reading 
RG6  6DW,  UK 

3  Rothamsted  Experimental  Station,  Harpenden,  Hertfordshire  AL5  2JQ,  UK. 


123 


DRAFT 


INTRODUCTION 

Computing  and  modelling  the  variogram  have  become  common  in  geostatistics, 
mainly  because  it,  or  its  equivalent  covariance  function,  is  a  prerequisite  for  local 
estimation  by  kriging.  It  has  merit  in  its  own  right,  however,  for  describing  spatial 
variation.  The  variogram  of  a  set  of  data,  usually  termed  an  experimental  or  sample 
variogram  to  distinguish  it  from  any  underlying  function,  summarizes  the  general 
nature  of  spatial  dependence  and  its  extent  and  also  the  magnitude  of  the  variation. 
Such  a  summary  may  then  be  interpreted  (cautiously)  in  terms  of  physical  process  - 
see,  for  example,  Webster  and  Cuanalo  (1975).  When  bounded  models,  such  as  the 
circular  and  spherical  functions,  fit  an  experimental  variogram  well  we  may  interpret 
them  as  representing  transition  structures;  they  suggest  patchy  distribution  of  values  in 
which  the  patches  have  extents  approximately  equal  to  the  ranges  of  the  models.  We 
can  take  this  a  step  further  and  look  for  boundaries  between  patches  at  this  scale,  as 
did  Webster  (1973,  1978)  and  Oliver  and  Webster  (1987).  We  know  that  in  computing 
the  variogram  we  lose  any  information  about  the  mean;  the  variogram  filters  it  out. 

In  the  last  few  years,  wavelets  have  started  to  find  application  in  the  earth  sciences 
(e.g.  Lark  and  Webster,  1999,  2001;  Oliver  et  al,  2000).  These  also  filter  data  and 
enable  structures  to  be  revealed  at  various  scales  by  multi-resolution  analysis. 

Wavelets  have  some  similarities  with  the  windowed  Fourier  functions,  but  the  wavelet 
transform  can  characterize  low  frequency  and  high  frequency  resolutions 
simultaneously.  Wavelets  operate  locally,  and  the  analysis  starts  with  the  choice  of  a 
mother  wavelet  w(x),  for  which  x  is  spatial  position,  that  can  be  dilated  or  shrunk. 

As  it  happens,  one  form  of  wavelet,  the  wavelet  of  Haar  (1910),  was  published  at 
about  the  same  time  that  Student  (in  Mercer  and  Hall,  1911)  realised  that  the  variation 
in  yield  between  adjacent  trial  plots  at  Rothamsted  Experimental  Station  was  derived 
from  two  sources:  an  autocorrelated  component  (structured)  and  a  random  component 
(locally  erratic).  This  idea  predates  the  variogram  of  geostatistics  by  several  decades, 
but  it  was  a  definite  precursor  to  it.  In  the  same  way,  the  Haar  wavelet  (Figure  1)  was 
described  long  before  recently  developed  wavelet  analysis  (Daubechies,  1988).  The 
Haar  wavelet  transform  is  remarkably  similar  to  the  experimental  variogram. 

Wavelets  can  filter  out  not  only  zero  and  first-order  moments,  but  also  higher-order 
moments.  They  can  also  separate  variation  into  components  from  several  spatial 
scales,  by  what  is  known  as  multi-resolution  analysis.  They  provide  a  much  richer 
repertoire  for  spatial  analysis  than  the  simple  variogram  because  they  are  not 
constrained  by  the  assumption  that  the  underlying  variable  is  random. 

The  purpose  of  this  paper  is  to  show  the  similarity  between  the  original  Haar  wavelet 
transform  and  the  variogram  of  a  set  of  data  and  then  to  suggest  an  extension  of  the 
latter  for  more  general  description. 


THE  WAVELET 

In  wavelet  analysis  a  signal  (information)  is  represented  in  terms  of  a  set  of  functions, 
which  we  denote  wsu(x),  where  5  signifies  scale  and  u  the  translation  space.  These  are 
basis  functions  or  kernels;  they  are  a  set  of  linearly  independent  functions  that  can  be 
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used  to  produce  all  admissible  functions  f(x)  (Strang  and  Nguyen,  1996).  Furthermore, 
they  are  constructed  from  the  single  mother  wavelet  w(x)  by  scaling  it  (by  2")  and 
translating  it  (by  u): 


W.yu(x)  = 


X 

U 


(1) 


Combining  the  basis  functions  with  the  wavelet  function  gives 


X  u 


(2) 


where  bsu  is  the  coefficient  of  the  basis  function  wyu(x)  at  scale  j,  translated  by  u,  and 
the  magnitude  of  which  provides  a  measure  of  the  significance  of  wsu(x). 

By  varying  5  and  u  we  can  analyse  a  signal  at  several  resolutions  in  a  multi-resolution 
analysis  of  the  wavelet  transform.  The  wavelet  function  w(x)  is  constructed  from  a 
scaling  function  <p(x)  which  can  be  specified  by  a  discrete  low-pass  filter  ln.  Similarly, 
the  wavelet  function  w(x)  can  be  specified  by  a  discrete  high-pass  filter  hn.  The  low- 
pass  filter  smoothes  the  signal,  while  the  high-pass  filter  retains  the  detail  in  the 
signal.  Daubechies  (1992)  and  Mallat  (1998)  describe  the  procedures  fully. 


THE  EXPERIMENTAL  VARIOGRAM 


We  start  with  the  familiar  formula  for  computing  the  variogram  based  on  the  method 
of  moments.  For  n  data  z(i),  i  =  1,2, . . . ,  n,  spaced  at  regular  intervals  along  a  line, 
Figure  la,  the  semivariance  is 


where  t,  an  integer,  is  the  lag  separating  the  data  points  z(i)  and  z(i  +  t).  By 


2  (n-t)  ,=I 


i  n-t 

=— y- 

n-tj^ 


z(i)  -  z(i  +  0 


incrementing  t  we  obtain  the  ordered  set,  the  experimental  variogram. 

To  see  the  relation  between  it  and  the  Haar  wavelet  transform  we  develop  the 
semivariance  in  terms  of  matrices.  Starting  with  t=  1,  for  which  the  comparisons  are 
shown  in  Figure  2a,  we  define  a  matrix  Bi  with  n  -  1  rows  and  n  columns: 
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1-10  0  0 
0  1  -10  0 


B,  =  —f= 


1 


0  0  0  .  .  .  0 


0 

0 


(4) 


Multiplying  B  i  by  the  data  vector  z=[z(  1 )  z(2)  . . .  z(n)]T  we  obtain 


'1-10  0  0 

.  .  o' 

"z(l)_ 

z(l)-z(2) 

0  1  -100 

.  .  0 

z(2) 

z(2)-z(3) 

0  0  0.. 

0  1  - 1 

_z(n)_ 

1 

=  42 

z(n-\)-  z(n)_ 

They'th  component  of  BjZ  is  {z(j)-z(j  +  1)}/V2.  So  to  obtain  7(1)  we  simply  square 
the  elements  of  Biz,  sum  them,  and  divide  by  n  -  1 : 


7(1)  =  “1t(B1z)t(B1z) 
n  —  1 


— - — zTB^B,z. 
n- 1 


(6) 


We  have  thereby  computed  the  square  of  the  Euclidean  norm  of  the  vector  Biz  and 
scaled  it  by  dividing  it  by  n  -  1.  The  square  of  the  Euclidean  norm  of  a  vector  z,  of 

length  n,  refers  to  the  following:  z1  z  =  ^”=]  (z,  )2 ,  which  is  the  same  as  the  sum  of  the 
square  of  the  components  of  z.  In  other  words,  Equation  (6)  expresses  f(  1)  as  a 
quadratic  form  -  see  also  Cressie  (1993),  page  73. 


Now,  to  compute  7(2)  we  define  a  matrix  B2: 


1  0  -1  0  0  ...  0 
01  0-10...  0 


1  0  -1 


(7) 


0  0 


0 
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Matrix  B2  has  n  -  2  rows  and  n  columns.  Multiplying  it  by  the  data  vector  z  we  obtain 


'10-100...  o' 
01  0-10...  0 

’zd)' 

z(2) 

z(l)-z(3) 

z(2)-z(4) 

1 

- 1 

•  1 

•  o 

•  o 

•  o 

•  o 

i 

_z(n)_ 

V2 

z(n-2)-z(ri)_ 

(8) 


As  for  f(l)  we  compute  f(2)  by 

7(2)  =  — ~(B2z)t(B2z)  =  ztBjB22  .  (9) 

n -2  n -2 

In  general  we  can  define  a  matrix  Bk  with  n-k  rows  and  n  columns  in  which  the  non¬ 
zero  values  1  and  —  1  are  placed  in  the  first  and  ( k  +  l)th  elements  of  the  first  row,  in 
the  second  and  (k  +  2)th  elements  of  the  second  row,  and  so  on.  On  multiplying 
through  we  obtain  the  general  expression  for  the  semivariance  equivalent  to  Equation 
(3)  as 


?(*)=■ 


-(Btz)T(Btz)  = 


n-k 


z‘B‘B,z. 


(10) 


THE  HAAR  WAVELET  TRANSFORM 


Let  us  now  look  at  the  original  Haar  wavelet  transform.  The  Haar  wavelet  transform 
splits  the  sequence  of  data,  usually  termed  the  signal  in  the  literature  on  wavelets,  into 
two  sets  of  components  corresponding  to  the  basis  coefficients  of  two  orthogonal 
subspaces.  One  component  contains  the  information  pertaining  to  the  local  averages 
of  the  signal,  the  smooth  or  low-frequency  components,  while  the  other  contains 
information  on  local  fluctuations  in  the  signal,  the  detail  or  high-frequency 
components.  The  Haar  wavelet  function  has  the  form  of  a  local  square  wave,  shown  in 
Figure  2. 
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We  proceed  as  follows.  We  suppose  that  the  signal  z  has  n  values  measured  at  equal 
intervals,  as  before,  but  with  the  constraint  that  n  is  a  power  of  2,  say  n=2m .  For 
reasons  that  will  become  clear  below  we  shall  designate  it  z\  =  z.  We  define  two 
matrices.  One,  which  we  denote  Li  containing  Is  and  Os,  is  for  the  low-frequency 
component  and  appears  in  Equation  (1 1).  The  matrix  Li  is  associated  with  the  scaling 
function  <p(x)  described  above.  The  other,  denoted  Hi  for  the  high-frequency 
component  and  containing  Is,  -  Is  and  Os,  occupies  an  analogous  position  in  Equation 
(12).  The  matrix  Hi  is  associated  with  wavelet  function  w(z)  described  above.  Both 
are  of  dimensions  n/2  rows  by  n  columns.  We  multiply  the  matrices,  scaled  by  1/*  2, 
by  the  data  vector  Z\  to  give 


=  LjZ,  = 


V2 


and 


d2  =  H,Zj  = 


S 


'1 

1  0 

0 

0 

0 

.  ...  0' 

'z,(l> 

Zi(l)  +  Zj(2) 

0 

0  1 

1 

0 

0 

.  .  .  0 

z,  (2) 

Zj  (3)  +  z,  (4) 

0 

0  0 

0 

1 

1 

.  .  .  0 

z,  (3) 

i 

Z]  (5)  +  Zj  (6) 

. 

. 

0 

0  0 

• 

• 

.011 

_Z\  in) 

- 

_Zi  (n  - 

- 1)  +  z,  (n)_ 

'1 

-1 

0 

0 

0 

0  .  . 

0  * 

~  z,  (1) 

Z,(l)-Z, 

0 

0 

1 

-1 

0 

0  .  . 

0 

z,(2) 

Zj (3) -  z 

0 

0 

0 

0 

1 

-1  .  . 

0 

z,(3) 

1 

Zj  (5)  —  z 

0 

0 

0 

.  .  0  1 

-lj 

_Z\  (n)_ 

Zj  (n  - 1)  - 

(ID 


.(12) 


We  thereby  obtain  a  vector  z2  containing  the  sums  of  neighbouring  values  (divided  by 
V2)  and  a  vector  d2  of  their  differences  (also  divided  by  V2).  The  elements  of  the  first 
provide  us  with  a  measure  of  the  smoothness  of  the  signal,  whilst  those  of  the  second 
tell  us  about  the  detail,  variation  or  change  in  Z|.  Note  that  z2  and  d2  consist  of  half  the 
number  of  components  of  z\,  which  is  due  to  the  number  of  rows  of  Li  and  Hi 
respectively. 

From  the  definition  of  Li  and  Hi  it  is  evident  that  the  rows  of  Lj  are  orthogonal  to  one 
another,  as  are  the  rows  of  Hi.  The  rows  of  Hi  are  also  orthogonal  to  those  of  L).  This 
is  true  for  all  orthogonal  wavelets.  Furthermore,  the  rows  of  Li  and  Hi,  respectively, 
form  an  orthonormal  set  of  basis  vectors.  From  z2  =  LjZi  and  d2  =  HjZi  (the 
coefficients  of  the  corresponding  basis  vectors)  we  can  compute  the  inverse  wavelet 
transform  to  obtain  zi  from  the  following 
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where  I  is  the  identity  matrix. 

We  can  now  take  the  analysis  a  step  further  by  splitting  z2  in  a  similar  way.  We  define 
two  matrices  L2  and  H2  (scaled  by  I/'s/ 2)  having  the  same  structure  as  Hi  and  Li,  but 
with  n/4  rows  by  n/2  columns  each  instead  of  nil  rows  by  n  columns,  and  we  multiply 
them  by  z2  to  give  the  vectors  z3=L2z2  and  d3=H2z2>  respectively.  As  in  Equation  (13), 
we  can  obtain  z2  from  z3  and  d3.  We  can  repeat  the  procedure  until  we  exhaust  the 
data.  At  each  step  we  split  the  smooth  component  of  the  signal  from  the  previous  step 
into  orthogonal  smooth  and  detail  components.  Finally  we  obtain  two  sets  of  vectors, 
Zj  and  dj,j=2,  3,...,  m  +  1,  where  n=2m.  This  series  of  transformations  will  provide  us 
with  a  means  of  analysing  the  signal  at  different  resolution  levels.  Furthermore, 
synthesizing  Zj  and  d j  for  levels  j= 2,  3,. . .,  m  +  1  to  obtain  zi  provides  for  a 
multiresolution  analysis;  see  Walker  (1999),  page  14. 

COMPARISON 

The  matrices  (of  the  successive  differences,  Equation  (5))  and  H,  (of  the  alternate 
differences,  Equation  (12))  show  the  similarities  between  the  variogram  and  the  Elaar 
wavelet  transform.  However,  there  is  not  a  one-to-one  correspondence  between  them; 
we  cannot  equate  the  variogram  at  any  particular  lag  with  any  high-frequency 
component  of  the  wavelet  transform.  This  is  because  in  computing  d j  at  any  particular 
j  we  have  lost  the  detail  in  the  original  signal  at  each  of  the  previous  j-1  steps,  and  in 
each  of  those  steps  we  have  halved  the  number  of  rows  in  the  matrices. 

Nevertheless,  by  expressing  the  variogram  in  terms  of  the  matrices  B*  in  Equation 
(10)  we  see  that  f(k)  does  contain  terms  corresponding  to  the  norm  of  the  high- 
frequency  components  H/Z;-  of  the  Haar  wavelet  transform. 

Consider  first  the  semivariance  y(l)  at  lag  t=  1  and  the  detail  d2=  H]Zi  of  the  wavelet 
transform.  Notice  how  matrix  Hi  is  related  to  matrix  Bj  used  to  form  the  y(l) —  it 
consists  of  the  odd-numbered  rows  of  Bi  only.  So  let  us  define  a  matrix  Ai  such  that 
Ai  =  Odd  rows  (Bi).  Then  for  the  case  t-\ 


(ZjH/HjZ,)  =  (z,A^A,Z|)  =  1  ^(z(2i~l)-z(2i) 

n- 1  ft  -1  ft- ly[  yfl 
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Following  from  the  above  we  find  that  the  semivariance  at  lag  t=  1  consists  of  two 
terms: 


1  n-l 

w>=^x  ^ 


z(i)  -  z(i  + 1) 


1  y?  T  z(2t  - 1)  -  z(2i) 
"n-lw  a/2 


+  ■ 


n/2-I 

7  X 

1  ;=i 


z(2/)-z(2i'  +  l) 


Odd  rows  of  B , 


Even  rows'  of  B , 


and  by  Equation  (14), 


=T_(Zxh,j,)+ 

n- 1 


1 


n  - 1 


n  /  2-1 

X 


;=i 


z(2z)-z(2i  +  l) 

V2  “ 


V 


In  words,  f(l)  contains  the  norm  of  the  detail  component  of  the  wavelet  transform  at 
level  1  scaled  by  l/(n-l).  When  t= 2  the  corresponding  f( 2)  contains  some 
components  of  the  norm  of  the  next  level  of  detail,  namely  d3=H2z2of  the  Haar 
wavelet  transform.  However,  the  norm  of  H2z2  contains  terms  that  are  not  in  f  (2) . 
The  same  is  true  for  other  values  of  t  and  different  levels  of  detail,  the  d,  of  the 
wavelet  transform. 

What  is  striking  about  the  matrices  B*  and  H,  is  that  both  consist  of  rows  the  elements 
of  which  are  only  0,  1  and  -  1.  In  the  Haar  wavelet  transform,  the  non-zero  scaling 
(smoothing)  coefficients  are  given  by  l\  =  /2  =1/V2,  whereas  the  non-zero  detail 
(difference)  filter  coefficients  are  h\  -  M'fl  and  h2  =  -  V'Jl.  All  the  other  coefficients 
are  zero:  /,  =  /r,  =  0,  for  all  i  >2.  The  wavelet  transform  is  obtained  by  shifting  by  2  the 
rows  of  the  matrices  L,  and  H,.  The  variogram,  however,  is  obtained  by  shifting  by  1 
its  rows  and  by  increasing  the  lag  between  the  non-zero  elements  of  each  its  rows. 
Figure  lb.  From  this  we  can  see  that  the  rows  of  B*  (a  running  difference),  which  are 
associated  with  a  two-point  variance,  are  related  to  filtering  data  with  the  Haar 
wavelet  filter  coefficients  1/V2  and  -1/V2.  In  general,  the  variance  of  n  data  z(i),  i  = 
1,2, . . . ,  n,  can  be  computed  by  a  circulant  matrix  the  rows  of  which  are  shifts  of  the 
following  elements:  ( n-\)/n ,  -\/n,. ..,  — 1  In.  Since  the  mean  of  z  is  given  by 

fl={z(l)+Z(2)+...+z(n)}fn,  and  the  variance  is  given  by  o2  - - l-LY  ■  We 

n  ~  1  i=i 

can  express  the  variance  as  (u  u)/(n-l)  which  is  the  same  as  (Pvz)  (Pvz)/(n-l), 
where 


(15) 
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'z(  1)-Jtx“ 

{(n  -  l)z(l)  -  z(2)  — ...  —  z(n)}/n 

z(  2)  -  /i 

_ 

{-  z(l)  +  (n~  l)z(2)  -  z( 3)  - ...  -  z(n)}/  n 

_z(n)~  H_ 

{-  z(l)  -  z( 2)  - ...  -  z(n  - 1)  +  in  - 1  )z(n)}/n 

P 


v 


(n-l) 

-1 


-1  -1 
(n-l)  -1 


-1 

-1 


-1  -1  ...  -1  (n-l) 


(16) 


It  is  not  difficult  to  show  that  (Pv)2  =  Pv  and  (PV)T=(PV).  In  other  words,  Pv  is  a 
projection  and  is  symmetric.  From  this  relation  we  can  rewrite  the  variance  of  z  as 
follows: 

62  =  (Pvz)T(Pvz)  _  (zTPvTPvz)  _  (zT(PvPv)z)  _  (zT(Pv)2z) 
n-l  n-l  n-l  n-l 


(zT(Pvz)) 
n  - 1 


This  means  that  the  variance  of  z  can  be  expressed  in  terms  of  a  projection  matrix  thus 
providing  a  very  simple  and  elegant  way  of  computing  such  a  measure.  Therefore,  we 
can  see  from  Equation  (17)  that  an  n-point  variance  can  be  couched  in  terms  of 
shifting  the  filter  coefficients  (n-l)/n,  -1/n,...,  -1/n. 


CONDITIONS  FOR  WAVELET  TRANSFORMS 

One  of  the  conditions  that  must  be  satisfied  for  the  wavelet  transform  is  that 


£/•,  =o, 


where  k  (finite)  corresponds  to  the  number  of  filter  coefficients. 
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Clearly,  from  the  paragraph  above,  this  holds  for  the  Haar  wavelet,  for  hi-  -h2  and  all 
the  other  coefficients  are  zero.  The  same  coefficients  are  also  present  in  the  formula 
for  the  variogram  but  in  a  different  order  because  of  the  lag  parameter  t  in  the 
variogram.  Nevertheless,  despite  their  differences,  both  formulae  provide  measures  of 
detail  or  variation  in  a  signal  at  several  scales. 


EXTENTION 

We  have  seen  above  that  there  are  just  two  non-zero  filter  elements  in  H,  and  L ,. 
Nevertheless,  wavelet  transform  filters  can  contain  more  elements.  So,  by  analogy, 
why  do  we  not  introduce  more  terms  into  the  variogram?  We  start  with 

f(t)  =  — —  Y.{h]z(i)  +  h2z(i  +  t)}2  ,  (18) 

n  —  t  ,=i 


where  h\  =  1  Nl  and  h2  =  -1/V2,  and  we  add  a  third  and  fourth  terms  to  give 

7(0  =  — - — YJ{h]z(i)  +  h2z(i  +  t)  +  h3z(i  +  2t)}2  ,  (19) 

and  n~2t  '=' 

7(0  =  • — - — V  {/I, z(0  +  h2z(i  +  0  +  h3z(i  +  2t)  +  h4z(i  +  3t)}2  ,  (20) 

n  -  3t  "f 


on  the  condition  that  at  least  we  have  =  0.  Note  that  in  Equation  (16),  the 

filter  coefficients  (n-\)/n,  -1  In  also  sum  to  0  (the  filter  has  one  vanishing 

moment).  Equation  (16)  shows  that  if  O1- 0  then  z=0  or  z=c  (constant). 

Furthermore,  in  general,  the  h\  and  hi  in  Equation  (19)  are  not  the  same  as  those  of 
either  Equation  (18)  or  Equation  (20).  Similarly,  the  h\,  h2  and  hi,  of  Equation  (20)  are 
not  related  to  those  of  Equation  (19),  and  so  on.  Increasing  the  number  of  filter 
coefficients  hj  allows  for  such  coefficients  to  satisfy  other  conditions.  We  shall 
elaborate  on  these  conditions  below. 

Remember  also  that  the  wavelet  transform  is  local.  And  so  we  should  expect  the  detail 
filter  coefficients  to  smooth  out  the  long-range  (low-frequency)  features  in  the  signal 
while  retaining  the  high-frequency  ones. 

Incorporating  these  coefficients  into  the  formula  for  the  variogram  may  be  beneficial 
in  that  it  embodies  more  information  than  just  the  differences  between  pairs  of  points. 
The  matrices  for  three  non-zero  coefficients  become 
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1 

J5- 

h2 

h , 

0 

0  .  . 

i 

O 

0 

to 

0  .  . 

.  0 

.  o 

_ i 

•  o 

•  o 

.  .  \ 

h2  - 

so  that 

7(  1)  ztB^B1z  , 

n  —  2 


(21) 


and 


\  0  h2  0  h3 

0  hx  0  h2  0  /z3 


0 

0 


0  0  0 


hv  0  h2  0  h3 


and 

7(2)  =  — 

n-4 


In  general  with  three  non-zero  coefficients  we  shall  have 

n-2k 


(22) 


This  is  more  restrictive  than  for  the  usual  variogram  because  we  have  only  n-2k 
terms  in  the  sum  instead  of  n-k,  and  so  we  run  out  of  data  more  readily. 
Alternatively,  we  may  see  it  as  demanding  more  data.  The  situation  will  be 
accentuated  if  we  want  more  than  three  filter  coefficients. 
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VANISHING  MOMENTS 


To  see  how  to  elaborate  the  variogram  as  suggested  above  we  need  to  understand  the 
concept  of  vanishing  moments.  We  explain  the  concept  via  the  y-transform  of  the 
filter  coefficients,  setting  y=exp(i  ft)),  where  i=V-l,  and  mis  the  frequency  variable, 
and  forcing  the  transform  to  have  77  zeros  at  co  =  0,  i.e.  y  =  1 .  (The  quantity  7]  is  the 
number  of  vanishing  moments,  and  27}  is  the  number  of  filter  coefficients).  This 
allows  the  filter  to  cancel  out  the  signal's  low-frequency  content.  We  shall  see  below 
that  cancelling  the  low  frequencies  is  related  to  cancelling  polynomials  of  order  up  to 
27]  -  1.  The  y-transform  of  the  filter  is  given  by 

277-1 

P^l^y-b 

k= 0 


or  equivalently 


277-1 

p2n  (<w)  =  £/i*exp(-ift*fc).  ('23') 

*=0 

(Note  that  p2v(y)  is  not  a  polynomial  in  y  since  the  exponents  of  y  are  negative).  For 
computational  convenience  the  filter  coefficients  are  given  by  ho,  h\,  ...,  h2T}- 1-  Forcing 
the  polynomial  equation  to  have  7]  zeros  at  co  =  0  means  that 

277-1 

p2rj(m)  =  ^hkexp(-icok)=  (O71  r(co),  (24) 

*= o 

where  r{(0)  is  a  polynomial  in  terms  of  ft).  Note  that  for  P2^(cd)  to  have  7]  zeros  at 
co  -  0  also  means  that  its  first  r\  -  1  derivatives  are  0  when  co  =  0.  That  is, 


=  XM-iO'expHfflO  =  0.  (25) 

d  C0m  *=o 

for  0  <  m  <  77  -  1 .  Factoring  (-  i)m  from  the  right  hand  side  of  Equation  (25)  and 
letting  m  =  0  we  obtain  the  vanishing  moment  conditions  for  values  of  m  where 
0  <  m  <  7]  -  1: 

277-1 

^ hkkm  =0,  which  is  the  same  as  (26) 

k=  0 


/70l  +  /j]l°  +h22°  +•••  +  ^27,-1  (2?7  -1)°  =  0,  when  m  =  0 
h0 O’  +  /ijl1  +h22]  +■■■  +  /72r)_,(277  -l)1  =  0,  when  m  =  1 

The 

fact  ho°n~l  +M’H  +h22n~l  +---  +  h2^l(2ri-\)T1~]  =0,  when  m  =  77 -1. 
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that  Equation  (26)  is  zero  means  that  increasing  the  number  of  vanishing  moments  in 
the  filter  will  also  reduce  to  zero  sections  of  length  2r\  of  polynomial  curves  of  degree 
0  <  m  <  rj  -  1 ,  when  the  filter  coefficients  are  convolved  with  such  polynomials.  That 
is,  not  only  will  the  filter  coefficients  /z,  reduce  to  zero  the  powers  of  0,  1,2,  . . 277—1, 
but  the  powers  of  (0+r),  (1+  r),  (2+  r ),  . . (2 77—  1  +  r),  where  r  is  any  fixed  real 
number,  are  also  reduced  to  zero.  It  is  important  to  note  that  Equation  (26)  is  not  a 
polynomial  function  of  k,  since  k  is  just  an  index  in  the  range  of  0<  k  <  277— 1. 

Interesting  properties  are  obtained  when  the  filter  coefficients  hi,  with  77  vanishing 
moments,  are  applied  to  a  polynomial  of  degree  m>rj.  This  operation  acts  on  the 
polynomial  like  a  derivative  of  order  m-77,  producing  a  polynomial  of  degree  m-r\. 

For  example,  consider  the  following  polynomial  s(y) 


s{y)=  (k  +  y)m, 

k= 0 


^zir1  _5L  f 


=  2*.X 

k= 0  i=0 


V1  ) 


y'k' 


(27) 


where  m  >  77,  hk  are  the  usual  filter  coefficients  with  exactly  77  vanishing  moments  and 
(, k+y)m  is  a  polynomial  in  y  of  degree  m.  Equation  (27)  can  be  rewritten  as 
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277-1 

where  ci  =  y£Jhkkm~l 

k= 0 


fO,  m-i<r]-~\<^i>m-r\  +  \ 

[non  zero  real  number,  m  -  i  >  rj  <=>  0  <  /  <  m  -  77. 


Clearly,  since  c,  is  nonzero  for  0  <  i  <  m-rj,  we  have  that  s(y)  =  X 


1=0 


c,-y 


is  a 


polynomial  in  y  of  degree  m-rj.  Therefore,  the  degree  of  the  polynomial  (k+y)m  has 
been  reduced  to  m-rj  from  m.  This  means  that  as  the  filter  mns  over  the  signal,  those 
portions  of  the  signal  consisting  of  polynomials  of  order  m  <  77  will  be  reduced  to  zero 
or  sustain  a  reduction  in  order  if  m  >  rj.  This  is  what  occurs  with  the  ordinary 
variogram  except  that  the  filter  coefficients  1/V2  and  —  1/V2  cancel  out  only  constant 
trend  (polynomials  of  order  zero).  Real  data  in  general  do  not  consist  exactly  of 
polynomial  functions  unless  a  polynomial  (exact  or  in  the  least  squares  sense)  is  fitted 
through  the  data  points.  Nevertheless,  as  the  filter’s  number  of  vanishing  moments 
increases,  the  wavelet-based  variograms  may  be  useful  to  identify  higher-order 
polynomial  global  trend  in  a  signal. 
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APPLICATION 

We  have  shown  above  the  relation  between  the  conventional  experimental  variogram 
and  the  filter  provided  by  the  Haar  wavelet  coefficients.  We  used  the  relation  to 
suggest  a  generalization  of  the  variogram  so  that  it  could  incorporate  other  filter 
coefficients  such  as  (but  not  limited  to)  Daubechies  s  2 77  coefficients.  Daubechies  s 
wavelet  filters  are  a  natural  extension  of  Haar’s  wavelet  filter.  For  an  even  number  of 
filter  coefficients  277,  the  wavelets  satisfy  77  orthogonal  conditions  as  well  as  77 
vanishing  moments. 

We  have  used  Daubechies's  high-frequency  filter  coefficients  for  277=  2,  4,  6,  and  8. 
(When  277  =  2  we  have  the  Haar  filter.)  Table  1  lists  their  values.  Since  the  wavelet 
transform  is  the  convolution  of  a  signal  with  a  wavelet  filter  followed  by  sub¬ 
sampling,  the  operation  requires  the  filter  (or  the  signal)  to  be  reversed.  The  values  in 
Table  1  are  listed  in  reverse  order  (relative  to  the  wavelet  transform)  to  show  the  order 
in  which  the  coefficients  have  been  applied  to  Equations  (18),  (29)  and  (30).  We 
emphasize  that  if  we  choose  a  different  order  for  the  wavelet  coefficients  we  shall 
obtain  different  results.  Nevertheless,  note  that  in  Equation  (3)  for  the  ordinary 
variogram  we  have  {h\z(i)+h2(i+t)}2={h2z(i)+  h\z(i+t)}2.  This  indicates  that  the  order 
is  not  important  since  there  are  only  two  filter  coefficients  that  differ  only  by  their  sign 
and  the  whole  quantity  is  squared. 

We  now  apply  the  ideas  to  soil  data  from  two  transects,  both  of  which  have  been  re¬ 
analysed  recently  by  Lark  and  Webster  (1999)  using  wavelets. 


Caragabal  transect 

Caragabal  is  on  the  Bland  Plain  of  Eastern  Australia.  The  topography  is  flat,  but  the 
Plain  is  pocked  by  gilgai,  that  is  by  small  depressions  typically  from  a  10  to  40  cm 
deep,  2  -  3  m  across,  and  20  to  40  m  apart  and  which  occur  in  a  seemingly  regular 
pattern.  The  soil  is  dominantly  clay  but  with  a  sandier  surface  horizon  of  variable 
thickness,  alkaline  and  locally  saline,  especially  in  the  gilgais  themselves. 

The  soil  was  sampled  by  Webster  (1977)  using  a  7.5-cm  diameter  corer  at  365  points 
at  4-m  intervals  along  a  transect.  We  used  the  records  of  the  soil’s  electrical 
conductivity,  a  measure  of  its  salinity,  at  30  -  40  cm  for  our  application.  Table  2 
summarizes  the  data.  The  measured  values  are  strongly  skewed,  and  we  have 
transformed  them  to  their  common  logarithms  for  our  analysis.  Table  2  summarizes 
their  statistics,  and  Figure  3  shows  the  (erratic)  trace  of  the  data  as  the  solid  line. 

Variograms  were  computed  with  the  usual  formula,  Equations  (3)  and  (18),  and  then 
with  its  wavelet-based  extension 
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f(t)  =  — —  V  fazQ)  +  h2z(i  + 1)  +  h3z(i  +  20  +  hAz{i  +  3r)}2  (29) 

n-3t~i 

and 

y(0  =  — - —  Y  {/jjzCO  +  h2z(i  +  0  +  h3z(i  +  20  +  hAz(i  +  30  +  ft5z(*  +  40  +  +  ^0}  >  (30) 

n-5t 

corresponding  to  Daubechies's  high-frequency  filter  coefficients  4  and  6,  respectively. 

Figure  4  shows  the  result.  The  three  variograms  are  remarkably  similar  to  one  another, 
though  we  see  that  as  the  number  of  filter  coefficients  increases  so  does  the  amplitude 
of  the  periodic  component.  See  Table  1  for  a  list  of  several  of  Daubechies’s  wavelets. 

This  is  a  consequence  of  the  wavelet’s  sensitivity  to  the  variation  in  the  data. 

Wavelets  in  general  eliminate  smooth  variation  such  as  trend  functions  that  can  be 
represented  by  low-order  polynomials,  but  not  the  local  fluctuation.  However,  the 
variation  introduced  in  Equations  (29)  and  (30)  by  the  additional  terms  z(i  +  kt )  and 
the  filter  coefficients  is  minimum. 


Sandford  transect 

Sandford-St-Martin  lies  on  the  Jurassic  scarplands  in  Central  England.  The  soil  there 
has  formed  by  the  weathering  of  the  gently  dipping  limestone,  sandstone,  siltstone, 
and  shale  of  the  Jurassic  era,  and  in  Recent  alluvium.  Webster  and  Cuanalo  (1975) 
surveyed  a  transect  3.2-km  long  across  the  outcrops  by  digging  pits  at  10-m  intervals 
and  recording  the  soil  to  1  m  depth.  For  present  purposes  we  have  analysed  the  clay 
content  at  two  depths,  topsoil  (centred  at  8  cm)  and  subsoil  (centred  at  65  cm),  in  the 
profile.  Their  statistics  are  summarized  in  Table  3,  and  Figure  5a  and  b  shows  the  data 
for  the  topsoil  and  subsoil,  respectively,  plotted  against  position  on  the  transect.  As  for 
Caragabal,  the  solid  lines  join  the  observed  values  and  show  much  point-to-point 
variation.  The  bold  dashed  lines  show  the  step  changes  in  the  mean  at  the  boundaries 
between  the  geological  formations. 

Again,  as  for  Caragabal,  we  computed  variograms  with  the  usual  formula  and  then 
with  its  extension  containing  the  wavelet-based  coefficients.  Figures  6  and  7  display 
the  results  for  the  topsoil  and  subsoil,  respectively.  The  variogram  of  the  clay  in  the 
topsoil  (Figure  6a)  seems  clearly  bounded.  Webster  and  Oliver  (1990)  modelled  it 
successfully  with  spherical  and  pentaspherical  functions,  with  a  sill  variance  of  about 
200  (%)2  and  a  range  of  about  250  m.  The  wavelet-based  variograms  with 
Daubechies's  high-frequency  filter  coefficients  4  and  6,  however,  increase  more 
gradually  and  have  smaller  maxima  -125  (%)2,  which  they  reach  at  shorter  lags 
(150  m  for  Daubechies’s  4  and  about  100  m  for  Daubechies’s  6,  Figure  6a).  All 
appear  to  approach  the  same  value  on  the  ordinate. 

To  explore  possible  causes  for  the  difference  between  the  usual  variogram  and  the 
wavelet  based  variograms  we  examined  the  influence  of  the  marked  boundaries  along 
the  transect.  We  removed  their  effect  by  computing  the  residuals  from  the  geological 
segment  means  and  computed  the  variograms  on  these  residuals.  The  result  is  the 
pooled  within-segment  variogram,  Figure  6b.  The  three  variograms  are  very  similar  to 
one  another.  The  usual  variogram  of  the  residuals  now  has  a  much  smaller  maximum 
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of  100(%)2  which  it  reaches  at  75  m.  This  is  the  same  for  the  wavelet-based 
variograms.  It  seems  that  the  latter  remove  the  effect  of  the  boundaries  (see  Figure  6a) 

_ Webster  and  Cuanalo  (1975)  identified  14.  Geostatisticians  might  regard  these 

local  step  functions  as  deterministic  features  to  be  removed  before  computing 
variograms  on  the  residuals.  It  appears  that  the  wavelet  variogram  filters  out  such 
trends  where  they  occur. 

The  variogram  of  the  subsoil  clay  content  was  fitted  best  by  a  pentaspherical  model.  It 
increases  to  a  sill  of  about  700(%)2  at  260  m  (the  range),  Figure  7a.  The  variogram 
from  Daubechies's  4  high-frequency  filter  coefficients  is  similar,  Figure  7a.  That  for 
Daubechies’s  6  rises  more  steeply  to  a  slightly  larger  maximum  at  about  the  same  lag 
distance  as  above,  Figure  7a.  The  variograms  computed  from  the  residuals  of  the 
segment  means.  Figure  7b,  are  similar  to  one  another. 

To  try  to  explain  the  reasons  for  the  difference  between  the  wavelet-based  variograms 
of  the  top-  and  sub-soil  we  did  an  analysis  of  variance  to  determine  the  proportion  of 
the  total  variance  accounted  for  by  geology.  The  results  of  the  analysis  of  variance  are 
given  in  Table  4.  The  proportion  of  the  topsoil  variance  explained  by  geology  is 
36.8%  and  that  for  the  subsoil  is  39%.  This  difference  is  not  sufficient  to  explain  why 
the  wavelet-based  variograms  filter  out  the  effects  of  the  boundaries  for  the  topsoil, 
but  not  the  subsoil.  A  possible  explanation  for  this  is  the  large  difference  between  the 
mean  for  the  first  part  of  the  transect,  positions  1  to  158,  and  that  of  the  second  part, 
positions  159  to  321. 


Simulated  transects 

To  attempt  to  understand  the  results  of  the  above  analyses  we  simulated  a  one¬ 
dimensional  sequence  of  values  in  three  ways  using  LU  decomposition  (see 
Goovaerts,  1997,  for  a  description  of  the  method).  We  used  a  spherical  variogram 
function  with  zero  nugget  and  a  sill  variance  of  1  for  each  simulation.  The  spherical 
function  is  widely  encountered  in  the  earth  sciences  and  provided  the  best  fit  to  the 
experimental  variograms  for  the  Sandford  data.  The  first  simulation  was  of  a  single 
process  based  on  this  variogram  with  a  range  of  45  units.  Figure  8a  (fine  line).  To 
emulate  the  effect  of  boundaries,  such  as  those  along  the  Sandford  transect,  the  second 
simulation  incorporated  a  stepped  function  superimposed  on  it  with  a  variance  of 
0.75,  Figure  8b.  The  lengths  of  the  segments  were  drawn  at  random  from  an 
exponential  process  with  a  mean  of  0.035.  Burgess  and  Webster  (1984)  discovered 
that  such  lengths  in  the  real  world  are  well  fitted  by  exponential  distributions. 

Although  we  simulated  the  steps  by  random  processes,  we  often  know  from  prior 
information  in  the  field  where  to  expect  them,  as  on  the  Sandford  transect.  They  can 
be  thought  as  a  form  of  regional  trend  because  they  result  in  predictable  changes  in  the 
mean. 

Regional  trend  also  results  from  a  gradual  change  in  values  with  autocorrelated 
random  fluctuation  superimposed.  An  example  of  this  is  the  top  of  the  sub  upper  chalk 
surface  described  by  Moffat  et  al.  (1986).  The  third  simulation  was  based  on  the  same 
variogram  function  as  the  first,  but  with  a  linear  trend  of  0.02  superimposed,  Figure 
8c. 
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Local  trend  or  drift  is  also  a  feature  of  the  real  world,  where  smooth  predictable 
change  in  a  property  occurs  over  short  distances.  The  electrical  resistivity  measured  at 
Beckesbourne  by  Webster  and  Burgess  (1980)  is  one  such  example.  Solifluction  and 
soil  creep  on  hillslopes  are  also  likely  to  give  rise  to  local  drift.  To  simulate  this  effect 
the  first  simulation  was  smoothed  by  an  1 1 -point  moving  average,  Figure  8a. 


Simulation  with  a  single  random  process 

The  sequence  of  simulated  values  in  Figure  8a  shows  them  fluctuating  about  a  mean 
value.  The  usual  variogram  and  the  wavelet-based  variograms  with  Daubechies’s 
high-frequency  filter  coefficients  4  and  6  were  computed  for  these  data.  Figure  9 
shows  the  result;  they  reach  a  similar  upper  bound,  but  the  wavelet  ones  fluctuate 
more  on  reaching  their  maximum. 


Simulation  with  a  stepped  function  superimposed  with  a  random  process 

Figure  8b  shows  the  simulation  with  stepped  functions  of  random  extent  (the  dashed 
line)  with  an  average  extent  of  45  units  and  variance  of  0.893,  superimposed  on  a 
random  process  generated  by  a  spherical  variogram  with  zero  nugget,  sill  of  0.25  and 
a  range  of  15  units.  The  latter  is  the  average  extent  of  the  within-segment  variation. 
There  were  nine  boundaries  or  step  changes  and  so  the  sequence  no  longer  fluctuates 
about  the  mean;  there  are  sections  where  the  average  value  is  small  and  others  where 
it  is  larger. 

The  usual  variogram  for  this  simulation  reaches  an  upper  bound  at  about  40  units, 
which  is  close  to  the  average  spacing  between  the  boundaries  of  45  units.  Figure  10a. 
There  is  a  kink  in  this  variogram  at  around  15  units  which  signals  the  nested  structure 
simulated.  The  wavelet-based  variogram  with  Daubechies's  high-frequency  filter 
coefficients  4  reaches  a  smaller  maximum  at  about  20  units  of  distance  and  then 
decreases,  as  does  that  for  Daubechies’s  6  which  fluctuates  more  about  the  upper 
bound.  These  results  suggest,  as  for  the  topsoil  clay  content  of  the  Sandford  transect, 
that  the  wavelet  variograms  filter  out  the  effects  of  the  long-range  transitions  in  the 
variation.  This  is  confirmed  by  the  variograms  computed  on  the  residuals  from  the 
class  means  for  this  simulation,  Figure  10b.  The  three  variograms  are  similar  to  each 
other  -  they  reach  an  upper  bound  of  about  0.25  which  is  reached  at  about  15  units. 
This  reflects  closely  the  short-range  function  used  for  simulation. 


Simulation  of  a  random  process  with  a  linear  trend  superimposed 

The  sequence  of  values  in  Figure  8c  shows  the  smooth  regional  progression  in  values 
from  let  to  right  with  random  fluctuation  superimposed  on  this.  Variograms  were 
computed  as  above  for  this  simulation.  There  is  a  clear  difference  between  the  usual 
variogram  and  the  wavelet-based  ones,  Figure  1 1 .  The  former  is  unbounded,  it  is  the 
kind  of  variogram  we  should  expect  where  the  data  contain  a  marked  linear  trend.  The 
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wavelet-based  variograms.  Figure  1 1,  describe  more  closely  the  spherical  function 
used  to  simulate  the  random  process.  This  result  also  suggests  that  the  wavelet 
variograms  filter  out  the  effects  of  regional  trend  in  the  data  as  well  as  those  of 
transition  features. 


Simulation  of  a  smooth  process 

Figure  8a  shows  the  1 1 -point  moving  average  of  the  underlying  random  process 
generated  from  a  spherical  model  with  zero  nugget,  a  sill  of  1  and  a  range  of  45  units 
(thick  solid  line).  The  local  smoothness  or  continuity  in  these  values  is  evident  in  the 
concave  slope  near  to  the  origin  of  the  usual  variogram,  Figure  12.  In  geostatistics 
such  data  should  be  analysed  using  a  structural  analysis  (Olea,  1975)  to  remove  the 
local  trend.  The  effect  of  the  latter  appears  to  be  amplified  in  the  wavelet-based 
variograms;  their  concave  slopes  near  to  the  origin  are  more  pronounced,  Figure  12. 
Thus  it  seems  that  wavelet-based  variograms  do  not  provide  a  solution  to  local  drift, 
and  this  should  continue  to  be  dealt  with  geostatistically. 


DISCUSSION  AND  CONCLUSION 

We  show  above  the  mathematical  similarity  of  the  method-of-moments  variogram  and 
those  derived  from  the  high  frequency  wavelet  coefficients,  Daubchies’s  4  and  6.  The 
experimental  variogram  computed  in  the  usual  way  is  sensitive  to  local  and  regional 
trends  in  the  data.  Their  presence  is  usually  evident  in  the  shape  of  the  variogram  and 
is  a  sign  that  the  stationarity  assumptions  on  which  much  of  the  standard  geostatistical 
technique  is  based  will  not  hold.  As  a  consequence  the  residuals  from  the  regional 
(arising  from  gradual  and  stepped  changes)  and  local  trend  should  be  determined  and 
used  for  further  analyses.  The  wavelet-based  variograms  appear  to  filter  out  regional 
trends  that  can  result  from  smooth  or  stepped  changes  in  the  regional  mean.  This 
would  avoid  the  need  to  examine  the  data  for  trend  first.  The  effects  of  local  trend  or 
drift,  which  are  less  common  in  most  earth  science  data,  appear  to  be  amplified  by  the 
wavelet-based  variograms,  however. 

This  is  a  preliminary  comparison  of  the  wavelet-based  variograms  with  Matheron’s 
estimator,  and  we  urge  caution  in  assuming  that  these  results  are  general.  More 
sequences  of  data  need  to  be  analysed  to  generalize  our  observations. 

Finally  we  should  add  that  our  circumspection  is  augmented  by  the  large  differences 
in  simulations  from  the  same  function,  but  with  different  seeds.  The  inconsistencies 
concern  us,  yet  many  researchers  continue  to  place  great  faith  in  simulation.  We  had 
some  difficulty  in  obtaining  sequences  with  experimental  variograms  like  those  with 
which  they  had  been  simulated. 
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Table  1.  High-frequency  Haar  and  Daubechies  filter  coefficients  (in  reverse  order) 
for  ?7=1,  2,  3,  and  4,  i.e.2?7=2  (Haar),  4,  6,  and  8. 


277  hi  h2 

h 

/l4 

h5 

h()  hj  h 

2  -0.7071  0.7071 

0 

0 

0 

0  0  0 

4  0.1294  0.2241 

-0.8365 

0.4830 

0 

0  0  0 

6  -0.0352  -0.0854 

0.1350 

0.4599 

-0.8069 

0.3327  0  0 

8  0.0106  0.0329 

-0.0308 

-0.1870 

0.0280 

0.6309  -0.7148  0.2304 

Table  2.  Summary  statistics  of  the  soil’s  electrical  conductivity  at  30-40  cm  at 

Caragabal  from  365  data 

Electrical  conductivity 

mS  cm'1 

Logio(mS  cm'1) 

Minimum 

0.06 

-1.214 

Maximum 

5.10 

0.707 

Mean 

0.985 

-0.230 

Median 

0.54 

-0.267 

Variance 

0.959 

0.192 

Standard  deviation 

0.975 

0.428 

Skewness 

1.64 

0.100 

Table  3.  Summary  statistics  of  the  clay  contenl 

l  as  percentage  by  weight  in  the  topsoil 

(5-10  cm)  and  subsoil  (62-67  cm)  at  Sandford  from  321  sampling  points. 

Topsoil 

Subsoil 

Minimum 

5.0 

0 

Maximum 

80.0 

90.0 

Mean 

25.6 

39.1 

Median 

20.0 

36.0 
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Variance 

Standard  deviation 
Skewness 


255.5 

16.0 

1.24 


936.8 

30.6 

0.16 


Table  4.  The  results  of  an  analysis  of  variance  for  the  top-  and  sub-soil  clay  content 
along  the  Sandford  transect. 


Topsoil  clay  content  (%) 

Source  of  variation 

Degrees  of 

Sum  of 

Mean 

Component 

freedom 

squares 

square 

of  variance 

Variation  due  to  geology 

14 

53013.4 

3786.7 

176.1 

Residual  variance 

306 

28744.5 

93.9 

93.9 

Total  variance 

320 

81757.8 

255.5 

255.5 

Subsoil  clay  content  (%) 

Variation  due  to  geology 

14 

187905.8 

13421.8 

622.7 

Residual  variance 

"306 

1 1 1 874 6 

_ 365.6 - 

320 

299780.4 

936.8 

936.8 

Total  variance 
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FIGURE  CAPTIONS 

Figure  1  (a)  Comparisons  made  in  computing  the  variogram,  (b)  Comparisons  for  the  Haar 
wavelet,  and  (c)  Points  embraced  when  computing  the  wavelet-based  variograms 

with  four  Daubechies’s  coefficients. 

Figure  2.  The  Haar  wavelet. 

Figure  3.  Trace  of  the  electrical  conductivity  (mS  cm'1)  transformed  to  common  logarithms 
at  Caragabal,  with  the  solid  line  joining  the  data. 

Figure  4  Variograms  of  the  electrical  conductivity.  The  usual  variogram  is  shown  by  •  , 
the  wavelet-based  variogram  with  Daubechies’s  4  coefficients  by  x,  and  with 
Daubechies’s  6  coefficients  by  □. 

Figure  5  Sandford  transect;  (a)  clay  content  in  topsoil  and  (b)  subsoil,  with  data  shown  by 

the  solid  lines  and  the  geologic  boundaries.by  the  dashed  bold  lines. 

Figure  6  Variograms  of  the  clay  content  in  the  topsoil  at  Sandford:  (a)  for  the  raw  data,  and 
(b)  for  the  residuals  from  the  means  of  the  geologic  classes.  The  usual  variogram  is 
shown  by  •  ,  the  wavelet-based  variogram  with  Daubechies’s  4  coefficients  by  x,  and 
with  Daubechies’s  6  coefficients  by  □. 

Figure  7  Variograms  of  the  clay  content  in  the  subsoil  at  Sandford:  (a)  for  the  raw  data,  and 
(b)  for  the  residuals  from  the  means  of  the  geologic  classes.  The  usual  variogram  is 
shown  by  •  ,  the  wavelet-based  variogram  with  Daubechies’s  4  coefficients  by  x,  and 
with  Daubechies’s  6  coefficients  by  □. 

Figure  8  Plots  of  the  simulated  data  using  LU  decomposition:  (a)  by  a  spherical  function 
with  zero  nugget,  unit  sill  and  a  range  of  45  units  (thin  line),  superimposed  by  a 
35-point  moving  average  shown  by  the  thicker  line,  (b)  by  a  spherical  function  with 
zero  nugget,  unit  sill,  range  of  15  units  and  a  random  stepped  process  with  an  average 
extent  of  45  units,  shown  as  the  heavier  dashed  line,  and  (c)  by  a  spherical  function 
with  zero  nugget,  unit  sill  and  a  range  of  45  units  with  a  linear  trend  of  0.02  added. 

Figure  9  Variograms  of  the  values  simulated  by  a  spherical  function  with  zero  nugget,  unit 
sill  and  a  range  of  45  units.  The  usual  variogram  is  shown  by  •  ,  the  wavelet-based 
variogram  with  Daubechies’s  4  coefficients  by  x,  and  with  Daubechies’s  6 
coefficients  by  □. 

Figure  10  Variograms  of  the  values  simulated  by  a  spherical  function  with  zero  nugget,  unit 
sill,  a  range  of  15  units,  and  a  step  function  with  an  average  extent  of  45  units:  (a)  for 
the  raw  data,  and  (b)  on  the  residuals  from  the  means  of  the  classes  created  by  the 
step  function.  The  usual  variogram  is  shown  by  •  ,  the  wavelet-based  variogram  with 
Daubechies’s  4  coefficients  by  x,  and  with  Daubechies’s  6  coefficients  by  0. 

Figure  11  Variograms  of  the  values  simulated  by  a  spherical  function  with  zero  nugget,  unit 
sill  and  a  range  of  45  units  smoothed  by  a  35-point  moving  average.  The  usual 
variogram  is  shown  by  •  ,  the  wavelet-based  variogram  with  Daubechies’s  4 
coefficients  by  x,  and  with  Daubechies’s  6  coefficients  by  0. 

Figure  12  Variograms  of  the  values  simulated  by  a  spherical  function  with  zero  nugget,  unit 
sill  and  a  range  of  45  units  with  an  added  linear  trend  of  0.02.  The  usual  variogram  is 
shown  by  •  ,  the  wavelet-based  variogram  with  Daubechies’s  4  coefficients  by  x, 
and  with  Daubechies’s  6  coefficients  by  □. 
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Fig  10 


